Article Text

Original research
Multimodal molecular landscape of response to Y90-resin microsphere radioembolization followed by nivolumab for advanced hepatocellular carcinoma
  1. Neslihan Arife Kaya1,
  2. David Tai2,3,4,
  3. Xinru Lim3,
  4. Jia Qi Lim1,
  5. Mai Chan Lau5,6,
  6. Denise Goh3,
  7. Cheryl Zi Jin Phua1,
  8. Felicia Yu Ting Wee3,
  9. Craig Ryan Joseph3,
  10. Jeffrey Chun Tatt Lim3,
  11. Zhen Wei Neo3,
  12. Jiangfeng Ye3,
  13. Lawrence Cheung3,
  14. Joycelyn Lee2,3,4,
  15. Kelvin S H Loke4,7,
  16. Apoorva Gogna8,
  17. Fei Yao1,
  18. May Yin Lee1,
  19. Timothy Wai Ho Shuen2,
  20. Han Chong Toh2,
  21. Axel Hilmer9,
  22. Yun Shen Chan10,
  23. Tony Kiat-Hon Lim11,
  24. Wai Leong Tam1,12,13,14,
  25. Su Pin Choo2,
  26. Joe Yeong3,5,11 and
  27. Weiwei Zhai1,15
  1. 1Genome Institute of Singapore (GIS), Agency for Science(A*STAR), Technology and Research, Singapore
  2. 2Division of Medical Oncology, National Cancer Centre Singapore, Singapore
  3. 3Institute of Molecular and Cell Biology (IMCB), Agency for Science, Technology and Research (A*STAR), Singapore
  4. 4Duke NUS Medical School, Singapore
  5. 5Singapore Immunology Network (SIgN), Agency for Science, Technology and Research (A*STAR), Singapore
  6. 6Bioinformatics Institute (BII), Agency of Science Technology and Research, Singapore
  7. 7Department of Nuclear Medicine and Molecular Imaging, Singapore General Hospital, Singapore
  8. 8Department of Vascular and Interventional Radiology, Singapore General Hospital, Singapore
  9. 9Institute of Pathology, Faculty of Medicine and University Hospital Cologne, University of Cologne, Koln, Cologne, Germany
  10. 10Guangzhou Laboratory, Guangzhou International Bio Island, Guangzhou, Guangdong Province, China
  11. 11Department of Anatomical Pathology, Singapore General Hospital, Singapore
  12. 12Cancer Science Institute of Singapore, National University of Singapore, Singapore
  13. 13NUS Centre for Cancer Research, Yong Loo Lin School of Medicine, National University of Singapore, Singapore
  14. 14Department of Biochemistry, Yong Loo Lin School of Medicine, National University of Singapore, Singapore
  15. 15Key Laboratory of Zoological Systematics and Evolution, Institute of Zoology, Chinese Academy of Sciences, Beijing, China
  1. Correspondence to Dr Weiwei Zhai; weiweizhai{at}; Dr Joe Yeong; yeongps{at}; Dr David Tai; david.tai.w.m{at}


Background Combination therapy with radioembolization (yttrium-90)-resin microspheres) followed by nivolumab has shown a promising response rate of 30.6% in a Phase II trial (CA209-678) for advanced hepatocellular carcinoma (HCC); however, the response mechanisms and relevant biomarkers remain unknown.

Methods By collecting both pretreatment and on-treatment samples, we performed multimodal profiling of tissue and blood samples and investigated molecular changes associated with favorable responses in 33 patients from the trial.

Results We found that higher tumor mutation burden, NCOR1 mutations and higher expression of interferon gamma pathways occurred more frequently in responders. Meanwhile, non-responders tended to be enriched for a novel Asian-specific transcriptomic subtype (Kaya_P2) with a high frequency of chromosome 16 deletions and upregulated cell cycle pathways. Strikingly, unlike other cancer types, we did not observe any association between T-cell populations and treatment response, but tumors from responders had a higher proportion of CXCL9+/CXCR3+ macrophages. Moreover, biomarkers discovered in previous immunotherapy trials were not predictive in the current cohort, suggesting a distinctive molecular landscape associated with differential responses to the combination therapy.

Conclusions This study unraveled extensive molecular changes underlying distinctive responses to the novel treatment and pinpointed new directions for harnessing combination therapy in patients with advanced HCC.

  • liver neoplasms
  • radiotherapy
  • immune checkpoint inhibitors

Data availability statement

Data are available in a public, open access repository. Raw whole exome sequencing (WES) and RNA sequencing fastq files are deposited in The European Genome-phenome Archive (EGA) under the study with accession code ‘EGAS00001006834’.

This is an open access article distributed in accordance with the Creative Commons Attribution Non Commercial (CC BY-NC 4.0) license, which permits others to distribute, remix, adapt, build upon this work non-commercially, and license their derivative works on different terms, provided the original work is properly cited, appropriate credit is given, any changes made indicated, and the use is non-commercial. See

Statistics from

Request Permissions

If you wish to reuse any or all of this article please use the link below which will take you to the Copyright Clearance Center’s RightsLink service. You will be able to get a quick price and instant permission to reuse the content in many different ways.


  • Immune checkpoint inhibitors (ICIs) have limited efficacy in treating advanced hepatocellular carcinomas (HCCs). Our previous phase II trial using yttrium-90 (Y90) radioembolization and ICI had shown a promising objective response rate, but molecular biomarkers for treatment response remain unknown.


  • For the first time, this study revealed multimodal biomarkers for Y90 radioembolization and anti-PD-1 (Programmed Cell Death Protein 1) combination treatment. These novel predictive biomarkers not only revealed new mechanisms responsible for the favorable treatment response, but also pointed to great potential for using these biomarkers harnessing the combination therapy for patients with advanced HCC.


  • Our study pointed new directions for using combination treatment in HCC with a unique set of biomarkers. Based on findings from this study, we can increase the response rate from ~15% in canonical ICIs to ~30% in this combination therapy. As many of the non-responders are enriched for an Asian enriched signature (Kaya_P2), we would expect patients from different ethnic backgrounds might have differential responses to this combination therapy, worth investigating in a future study.


Hepatocellular carcinoma (HCC) is the most common form of liver cancer and is the third leading cause of cancer-related deaths worldwide.1 Most patients with HCC are diagnosed at an advanced stage when systemic therapy is indicated.2 The landscape of systemic therapy in HCC has changed considerably over the past decade. It began with the use of sorafenib, a multitargeted tyrosine kinase inhibitor3 4 followed by others of a similar class such as lenvatinib,5 regorafenib,6 and cabozantinib.7 With advances in immunotherapy, monotherapy with immune checkpoint inhibitors (ICIs), such as nivolumab and pembrolizumab, showed early promise in terms of efficacy in single-arm studies; however, primary endpoints were not met in subsequent randomized trials.8 Combination therapy with atezolizumab plus bevacizumab showed superior efficacy compared with sorafenib (30% vs 11% objective response rate (ORR) and 19.2 vs 13.4 months median overall survival, respectively) and has become the new first-line treatment for advanced HCC.9 More recently, the combination of a single priming dose of tremelimumab and durvalumab displayed superior efficacy and fewer adverse effects compared with sorafenib and is thus positioned as an alternative option in first-line treatment for advanced HCC.10 Despite rapid progress in the treatment of HCC, response rates in advanced HCC are still modest and grade 3/4 treatment-related adverse effects remain a challenge, occurring in 20–40% of cases.11–13 Novel treatment strategies, personalized medicine and alleviation of treatment-related adverse events remain unmet needs in HCC.

Therapeutic synergy refers to an increase in the efficacy of a combination treatment compared with that achieved with any of its individual components alone.14 15 One of the obvious limitations of ICIs is that a large proportion of advanced HCCs are immune-excluded (ie, ‘cold’ tumors) and do not respond to checkpoint inhibition. Locoregional treatments include yttrium-90 (Y90) radioembolization (RE), which is a radiation therapy known to activate the immune system.16 Therapeutic synergism has been observed when combining ICIs and radiotherapy (RT) in preclinical HCC models, possibly through interferon gamma (IFN-γ)-STAT3 signaling pathway-mediated upregulation of PD-L1 (Programmed Death-Ligand 1).16 To evaluate the efficacy of combination therapy with ICI and RT in patients, a Phase II trial (CA209-678) was conducted to investigate the efficacy and safety of Y90-RE followed by nivolumab in patients with advanced HCC (NCT03033446, n=36).17 Interestingly, the ORR was found to be 30.6% in all patients and 43.5% when restricted to patients with only intrahepatic tumors, suggesting possible synergy between Y90-RE and ICI. Despite the encouraging treatment response, the underlying molecular mechanisms and biomarkers that can predict the clinical response remain unknown.

In this study, we collected longitudinal (pretreatment and on-treatment) biopsies and blood specimens obtained before Y90 treatment and after the first dose of nivolumab from patients included in the Phase II trial (CA209-678). By whole exome sequencing (WES), RNA-sequencing (RNA-seq), and multiplex immunohistochemistry (mIHC) of tissue biopsies as well as Luminex assays of patient blood samples from 33 patients, we systematically investigated molecular mechanisms associated with the favorable response to the combination treatment. As multimodal perusal of dynamic changes in patient matched pretreatment and on-treatment tumor biopsies remain rare in advanced HCC. Our study provided one of the largest data set until such efforts allowing study of mechanistic insights behind combination Y90-RE and nivolumab. By comparing genomic, transcriptomic and tumor microenvironmental differences between responders and non-responders, we discovered extensive molecular phenotypes delineating differential responses to the combination therapy and pinpointed new directions for harnessing combination therapy in patients with advanced HCC.


Sample collection and genomic sequencing

Thirty-six patients from the trial received a single dose of intrahepatic arterial Y-90 RE followed by intravenous nivolumab 21 (±3 days) after RE. Intravenous nivolumab (240 mg) was continued every 2 weeks thereafter until the occurrence of disease progression, death, severe toxicities or at the physician’s discretion. For each patient, pretreatment and on-treatment biopsies were collected from the same liver lesion before Y90-RE and after one dose of intravenous nivolumab (figure 1A). Based on the RECIST (Response Evaluation Criteria in Solid Tumours) V.1.1 guidelines, the ORR was found to be 30.6%, with 1 case of complete response (CR), 10 of partial response (PR), 11 of stable disease (SD), and 11 cases of progressive disease (PD). Three patients were non-evaluable. To compare treatment responses, patients who achieved either CR or PR were grouped as ‘responders (R, n=11) and patients who had either SD or PD were defined as ‘non-responders’ (NR, n=22) (figure 1A). To identify genomic changes associated with treatment response, we successfully performed WES (n=29), RNA-seq (n=28), and mIHC (n=29) of tumor biopsies obtained before and during treatment. In addition, the expression of a set of liquid biomarkers was determined by Luminex assays (n=33) of blood samples obtained before and during treatment (figure 1A, online supplemental figure 1).

Supplemental material

Figure 1

Study design and baseline genomic differences between responders (R) and non-responders (NR). In this figure, patients with baseline DNA data was used for comparisons (7 R, 15 NR, n=22). (A) Schematic diagram of the trial design (top) and summary of biospecimen collection distribution (bottom). Gray boxes indicate missing data where biopsy or blood sample collection was not successful or samples did not meet the quality control criteria. Number of patients with only baseline (pre) and paired samples are indicated on the right side of the heatmap for each assay. (B) Non-synonymous tumor mutation burden (TMB) comparison between R and NR (Methods). (C) Baseline clonal TMB comparison between R and NR. (D) Baseline subclonal TMB comparison between R and NR. (E) Total baseline SNP neoantigen count comparison between R and NR. (F) Landscape of driver mutations in the cohort. Baseline frequencies of driver genes in R and NR at the patient level are shown in barplots (bottom right). Genomic changes that are significantly different in the pretreatment samples between R and NR are shown in red. F, female; M, male; HBV, Hepatitis B Virus; HCV, Hepatitis C Virus; SNP, Single Nucleotide Polymorphism; CNV, Copy Number Variation; mIHC, multiplex immunohistochemistry; NeoAg, neoantigen; RNA-seq, RNA-sequencing; WES, whole exome sequencing; Y90, yttrium-90.

A collection of genomic changes predict response at the baseline

Among the 33 patients, we successfully obtained WES data for 65 baseline and on-treatment biopsies from 29 patients with a mean sequencing depth of 125×. After removing low purity samples (n=9) and one hypermutated sample, 27 patients with both pretreatment (n=22) or on-treatment (n=22) samples were available for further analysis. As TMB is one of the most important genomic predictors of ICI responses across many cancer types,18 but was not indicative in HCC,19 we compared the non-synonymous TMB between R and NR, surprisingly, a significantly higher TMB was found in R in the baseline samples (p=0.044, figure 1B). This difference remained significant after controlling for confounding variables such as tumor purity (online supplemental figure 2). In earlier studies, clonal TMB was found to be a reliable predictive biomarker of ICI responses in other cancer types.20 When we conducted a separate comparison of the clonal and subclonal TMB of pretreatment samples, we found that higher clonal TMB was a more significant predictor of patient response than total TMB (p=0.017, figure 1C). In contrast, subclonal TMB was not associated with response (p=0.34, figure 1D). When we compared mutations that are likely to produce neoantigens, we found that the number of neoantigens was significantly higher in R compared with NR (p=0.012, figure 1E). Although TMB has not been associated with response in previous immunotherapy trials in HCC,19 our findings indicate that it is a significant predictor of patient response to the combination therapy.

In addition to multiple factors related to TMB, we evaluated potential associations of oncogenic driver mutations with treatment response. By compiling several lists of cancer driver genes in HCC (Methods) and comparing their mutational frequencies between R and NR (figure 1F), we found that canonical driver genes such as CTNNB1 (36%) and TP53 (45%) were not associated with treatment response. Strikingly, NCOR1, which encodes nuclear receptor co-repressor 1 protein, was mutated in three out of the seven R, but was not observed in NR (0/15) (p=0.02, figure 1F, online supplemental figure 3). As a gene that mediates ligand-independent repression of thyroid hormone and retinoic acid receptor transcription, inactivating mutations in NCOR1 have been reported to be associated with better prognosis in prostate and bladder cancers.21 22 Interestingly, all three NCOR1 variants in R were truncating mutations (one frameshift and two stop codons) indicating that NCOR1 inactivation may be associated with patient responses to the combination therapy.

In addition to TMB and driver mutations, previous pan-cancer studies revealed a strong correlation between genome instability and immune evasion.23 Moreover, patients with higher somatic copy number alteration (SCNA) levels tend to be associated with poorer response in immunotherapy.23 However, when we compared the genomic alterations based on the Genomic Instability Index (GII) and SCNA score between R and NR, there were no significant differences in the genome instability of both chromosomal and focal alterations between two groups (online supplemental figure 4A–C). Interestingly, when we analyzed the overall genome instability and compared individual chromosome arms, we identified significant differences in chromosome 16p deletions between R and NR, with all the deletions found in NR (n=9/15, p=0.02, figure 1D, online supplemental figure 5). Chromosome 16 deletions were recently found to be associated with a new transcriptomic subtype (P2, denoted as Kaya_P2) enriched in Asians.24 The Kaya_P2 subtype was found to have high alpha-fetoprotein (AFP) values and is immunologically suppressive, possibly leading to extremely poor treatment response in NR (see later sections and Discussion).

Strong clonal turnover in R drives favorable treatment response

We observed that the response to the combination therapy was predicted by TMB and neoantigen burden as well as NCOR1 mutations, while NR was strongly associated with chromosome 16 deletions in baseline samples. In addition to baseline differences, genomic changes associated with treatment response also show good predictive ability.25 When we compared the clonal and subclonal TMB between pretreatment and on-treatment biopsies, we indeed found a significant decrease of clonal TMB in R, but not in NR following the treatment (figure 2A). Neoantigen burden was also significantly decreased in R (figure 2B). Thus, the combination treatment seemed to have triggered the loss of multiple mutations from the primary tumor.

Figure 2

Genomic changes during treatment determined by comparing pretreatment and on-treatment samples. (A) Changes in clonal TMB in responders (R, left) and non-responders (NR, right). (B) Changes in the number of neoantigens in R (left) and NR (right). (C) Landscape of different types of mutations (contracting, persistent or expanding) among R and NR using patients with paired pretreatment and on-treatment biopsies (7 R, 10 NR, n=17, see figure 1A). Mutational status of driver genes is also indicated below. (D) The proportion of contracting mutations in R and NR using patients with paired pre and on treatment biopsies (7 R, 10 NR, n=17, see figure 1A). (E) Fish plot showing clonal dynamics under treatment for a representative responder (P012). (F) Fish plot of a representative non-responder (P011). (G) The proportion of contracting neoantigens (NeoAg) in R and NR using patients with paired pre and on treatment biopsies (7 R, 10 NR, n=17, see figure 1A). On, on-treatment; Pre, pretreatment; TMB, tumor mutation burden.

To further understand the clonal dynamics during treatment, we analyzed the genomic evolution using paired pretreatment and on-treatment biopsies (7 R and 10 NR). Using an approach similar to a previous study by Riaz et al,25 we categorized each mutation based on cancer cell fraction (CCF) as ‘contracting’ (≥10% decrease in CCF, see Methods), ‘expanding’ (≥10% increase in CCF) and ‘persistent’ (<10% difference) based on changes in their CCFs. Indeed, we found a higher proportion of contracting mutations in R compared with NR (figure 2C–E). Conversely, persistent and expanding mutations were dominant in the majority of NR (figure 2C,F, online supplemental figure 6), suggesting very minor treatment effects in NR.

Given the marked clonal changes in R, we explored possible mechanisms driving favorable responses. As RT can lead to ‘leakage’ of neoantigens and subsequent elimination of antigen-positive clones by the immune system,25 we compared the percentage of contracting neoantigens in R and NR. In accordance with the decreased neoantigen burden in R, the percentage of contracting neoantigens was also significantly higher in R compared with NR (p=0.036, figure 2G). Taken together, the clonal dynamics revealed here indicate that immune activation in R might lead to ‘tree pruning’ dynamics, whereby multiple branches of the tumor phylogeny are sheared away during the treatment. On the other hand, NR tended to be dominated by relatively stable, or even expanding clonal trajectories.

Inflammation pathways and a novel transcriptomic subtype predict patient response

The genetic comparisons provided important insights into the genomic changes associated with treatment response in our cohort. We then investigated the corresponding phenotypic differences associated with treatment response. Transcriptomic comparison of baseline (pretreatment) samples revealed 157 upregulated and 126 downregulated genes in R compared with NR (figure 3A, online supplemental figure 7 and table 1). The upregulated genes were enriched in several pathways such as the IFN-γ pathway (IFNGR1) and epithelial-to-mesenchymal transition (RASGRF1 and COL23A1) (figure 3A,B). On the other hand, genes enriched in NR were more active in cell cycle-related pathways including MYC and E2F targets (figure 3A,B).

Figure 3

Baseline transcriptomic differences between responders (R) and non-responders (NR). In this figure, patients with baseline RNA-seq data was used for comparison (8 R, 16 NR, n=24, see figure 1A). (A) Volcano plot of differentially expressed genes (DEGs) between R and NR at the baseline. (B) Significantly enriched pathways based on DEGs between R and NR. Pathways with positive scores are significantly enriched in R and pathways with negative scores are enriched in NR. (C) Enrichment of known HCC transcriptomic subtypes in R and NR. (D) Comparison of subtype signature scores (Hoshida S1) between R and NR. (E) Comparison of subtype signature scores (Hoshida S2) between R and NR. (F) Comparison of subtype signature scores (Kaya P2) between R and NR. (G) Comparison of subtype signature scores (Yamashita EpCAM Up) between R and NR. (H) Significance of previously reported signatures of immune checkpoint blockade in R and NR (n=22). HCC, hepatocellular carcinoma; RNA-seq, RNA-sequencing.

In earlier studies of HCCs, several different transcriptomic subtypes were identified in patients from different ethnic backgrounds.24 26–28 An Asian enriched subtype (P2, denoted Kaya_P2) linked to poor overall survival was found to be associated with upregulation of proliferation pathways, such as MYC targets, and chromosome 16 deletion.24 To understand the possible association between patients with different transcriptomic subtypes and treatment response, we used gene signatures from known HCC subtypes and calculated enrichment scores across samples (Methods, online supplemental table 2). The highly enriched pathways identified in NR included signatures of proliferative HCC subtypes such as Kaya_P2,24 Hoshida S227 and Yamashita’s EpCAM signature29 (figure 3C–G).

Supplemental material

In addition to the biomarkers identified in this study, the increased application of immunotherapy has led to the identification of a large collection of biomarkers that can be used to stratify patient response.30–32 Therefore, we curated a set of 22 signatures from the literature, including markers reported for a cohort receiving immunotherapy for HCC,32 TIDE (Tumor Immune Dysfunction and Exclusion),30 radiotherapy response signatures33 34 and a nine-gene WNT (Wingless/integrated) pathway score.31 When we compared the expression levels of these biomarkers in R and NR at the baseline, we found that none were significant for stratification of patients in this cohort (figure 3H). Furthermore, we did not observe any differences when comparing these biomarkers separately for non-viral and HBV+ patients as well as patients without multiple disease (extrahepatic spread or multiple liver tumors) (online supplemental figures 8–9). Thus, the molecular phenotypes associated with combination treatment might be different from those associated with single-agent immunotherapy.

Differential transcriptomic changes delineating the treatment response

At the baseline, we observed that IFN-γ pathways tended to be highly expressed in R, while cell cycle related pathways were enriched in NR. To further dissect the phenotypic changes associated with treatment response, we compared the transcriptomic profiles of baseline and on-treatment samples. Comparison of the fold-changes in gene transcription during the treatment response showed similar trends in the expression of most genes in R and NR, with the majority of differentially expressed genes (DEGs) (92%) changing in the same direction (figure 4A). Interestingly, only 8% of DEGs showed changes in different directions (up in R, but down in NR or vice versa, figure 4A). When we identified DEGs by comparison of pretreatment and on-treatment samples in R and NR separately, we observed both convergent and divergent responses in R and NR (figure 4B). For example, immune-related pathways, such as IFN-γ response, were upregulated in both R and NR. Similarly, cell cycle pathways, such as MYC targets, were downregulated after the treatment in both groups. Interestingly, several important metabolism-related pathways, such as bile acid metabolism and peroxisomes, were downregulated in R, but not in NR (figure 4C, online supplemental figures 10 and 11), indicating potential divergence in the mechanism of the treatment response in the two groups.

Figure 4

Convergent and divergent treatment responses between responders (R) and non-responders (NR). (A) Log fold-change in differentially expressed genes (DEGs) between on-treatment and pretreatment samples in R and NR. (B) Summary of the number of DEGs between pretreatment and on-treatment samples in R and NR. (C) Pathways that are significantly upregulated and downregulated with treatment in both R and NR. Stars indicate a significant enrichment (fdr<0.05). While the majority of significant pathways showed similar direction in both R and NR, metabolism related pathways such as bile acid metabolism were changed in opposite directions in R and NR. Pathways changing in different directions are shown below the horizontal dashed line. (D) Significance of the comparison of reported immune checkpoint blockade response signatures (n=22) before and after treatment in R and NR. (E) Boxplot of selected signatures. Changes in expression are shown in R and NR separately. HCC, hepatocellular carcinoma; fdr, false discovery rate; On, on-treatment; Pre, pretreatment.

To compare the immune-related changes in more detail, we determined the differences in immunotherapy response signatures (figure 4D)30–34 and observed significant changes in several signatures that occurred only in R (figure 4E, online supplemental figure 12). Thus, although immune response pathways changed in similar directions in both R and NR, a greater increase in immune response signatures (eg, cytolytic activity), was often observed in R. In contrast, previous studies showed that the decrease in WNT signaling is often associated with increased IFN-γ signaling by T cells.31 35 Indeed, we found a significant decrease in the nine-gene WNT score in R, but not in NR (figure 4E). These findings imply that greater changes in immune-related pathways (ie, increase) and the WNT pathway (ie, decrease) are associated with a better treatment response in R.

Macrophages, but not T cells, are associated with the favorable response

From the genotypic and phenotypic analyses, we gained a good overview of the genomic and transcriptomic changes associated with the treatment responses in this cohort. Since the tumor microenvironment (TME) can strongly affect treatment response, we further explored its composition by immune deconvolution of the bulk RNA as well as mIHC of the tumor tissues. Using a well-established immune deconvolution method,36 we estimated the proportion of 14 immune cell types from the bulk expression data and clustered tumor samples into ‘hot’ and ‘cold’ subgroups (see Methods). Although we observed a higher infiltration of immune cells in the post-treatment sample (figure 5A), unlike other cancer types, we did not observe any significant enrichment in immune ‘hot’ tumors or individual lymphocyte cell scores, including CD8+ T cells, in R at the baseline (figure 5A). However, when we compared the non-lymphocyte cells between R and NR at the baseline (online supplemental figure 13), we observed significant enrichment in macrophages (p=0.02) and mast cells (p=0.03) in R (figure 5A).

Figure 5

Comparison of the tumor microenvironment (TME) in responders (R) and non-responders (NR). (A) Heatmap of cellular scores derived from bulk RNA-seq deconvolution using the method described by Danaher et al.36 Bar plots on the right show effect sizes when comparing R and NR using pretreatment (left) and on-treatment samples (right) separately. Red bars indicate a higher score in R and blue bars indicate a higher score in NR. Borders around the bars indicate significant p values. (B) Comparisons of the percentages of different cell types (identified by cellular markers) identified in the mIHC assay. Circle color represents the effect size when comparing marker values between R and NR. Black borders around the circles indicate significant p values. (C) Boxplots of markers showing significantly different expression between R and NR at the baseline using patients with available baseline mIHC data (7 R, 20 NR, n=27, see figure 1A). (D) mIHC images of the same significantly different markers in representative samples. Square boxes indicate positive cells for the marker at the header of each panel. IHC, immunohistochemistry; fdr, false discovery rate; mIHC, multiplex IHC; On, on-treatment; Pre, pretreatment; RNA-seq, RNA-sequencing; TIL, Tumor-Infiltrating Lymphocyte; NK, Natural Killer; Th1, T helper type 1; DC, Dendritic Cell, ES, Enrichment Score.

To further explore the role of the TME in treatment response, we conducted mIHC with a panel of immune lineage markers of T cells (eg, CD8) and macrophages (eg, CD68), as well as immune checkpoint or activation markers (PD-L1, PD-1, LAG-3 and CD38). Markers (CXCL9-CXCR3) reported to have good predictive response power in multiple solid tumors treated with immunotherapy20 37 38 were also included in the panel, although no differences in the expression of immune checkpoint molecules, such as PD-1, LAG-3 and PD-L1, were observed at the baseline (figure 5B). In accordance with the bulk deconvolution results, there were no differences in the baseline percentage of CD8+ T cells between R and NR (figure 5B), while there was a higher baseline percentage of CD68+ cells in R (figure 5B–D). Intriguingly, CD68+ macrophages and CD45+ immune cells, which also co-express CXCL9 and the receptor CXCR3,20 39 demonstrated good predictive power for differentiating R from NR (figure 5C,D). Further analysis showed that CXCL9+ cells in the TME were predominantly CD68+ macrophages (online supplemental figure 14).

In addition to characterization of the TME in the tumor tissue, plasma cytokine analysis can provide additional insights into the physiological and pathological processes that could be used to aid diagnosis and treatment.40 By leveraging the Luminex technology which enables simultaneous analysis of 65 cytokine targets, we compared the extent of systemic inflammation in pretreatment and on-treatment blood samples of the patients. Among the pretreatment samples, CXCL9, eotaxin-3 and LIF were found to provide predictive value for distinguishing R from NR (figure 6A). LIF is a pleiotropic cytokine that generates a local immunosuppressive microenvironment through various mechanisms, such as by promoting radioresistance41 and chemoresistance.42 Elevated levels of circulating LIF correlate with tumor recurrence and contribute to chemoresistance.42 Tumors with high levels of LIF tend to have tumor-associated macrophage (TAM) infiltration with repressed CXCL9 induction as well as CD8+ T-cell infiltration.43 Eotaxin-3 has also been associated with tumorigenesis and progression of HCC through the recruitment of immunosuppressive TAMs and myeloid-derived suppressor cells.44–46 These reports lend support to our observation of reduced CXCL9 with increased eotaxin-3 and LIF levels in the blood, as well as the higher baseline percentage of CD68+ cells in the TME of R pretreatment (figure 5B–D). Among the on-treatment samples, several immune-related markers, such as CXCL11, IL-12p70, IL-1β and CCL2, showed predictive value for distinguishing R from NR (figure 6B). Overall, these findings indicate that patterns of immune-related markers in the TME and multiple cytokines in the blood can provide important predictive value for distinguishing R from NR (figure 6C).

Figure 6

Comparison of blood biomarkers between responders (R) and non-responders (NR). Markers showing a significant difference in expression at either (A) pretreatment or (B) on-treatment between R and NR are shown. All 33 patients have available data for Luminex assay (see figure 1A). (C) Summary of baseline biomarkers identified in our study across all patients (n=33). Progression-free survival (PFS) times of patients are annotated as bars at the top. Response and extrahepatic spread status of patients are also annotated at the top. mIHC, multiplex immunohistochemistry; TMB, tumor mutation burden.


Using patient samples from the Phase II trial (CA209-678),17 we systematically investigated molecular changes associated with differential response to treatment with Y90-RE followed by nivolumab. Through WES, RNA-sequencing, and mIHC of tumor biopsies as well as Luminex assays of blood samples obtained before and during treatment, we identified higher TMB and neoantigen burden as well as NCOR1 mutations as genetic changes enriched in R at the baseline. Phenotypically, tumors from R showed more inflammation with enrichment of IFN-γ pathways and a greater proportion of macrophages. In contrast, NR tended to be enriched for a novel Asian-specific transcriptomic subtype (ie, Kaya_P2) with a high frequency of chromosome 16 deletions and upregulated cell cycle pathways. Comparison of the genomic differences between baseline and on-treatment biopsies showed that R tended to have greater clonal changes manifested as a significant decrease in the TMB and neoantigen burden at the genetic level. Despite a convergent upregulation of immune response and downregulation of cell cycle pathways in both R and NR, R had a more pronounced increase in several immune signatures, including IFN-γ and cytolytic activity, but a significant decrease in WNT signaling. Luminex assays of the blood samples revealed strong segregation of several chemokines, including CXCL9, CXCL11 and IL-1β, between R and NR. Taken together, we have unraveled extensive molecular changes that delineate differential responses to combination therapy and pinpointed new opportunities for harnessing combination therapy in advanced HCC.

Our phenotypic analysis of pretreatment samples showed that there were no differences in the proportion of hot and cold tumors between NR and R, although the IFN-γ (IFNGR1) pathway was upregulated in R with a higher percentage of CXCL9+/CXCR3+ CD68+ macrophages and CD45+ immune cells. These findings highlighted the predictive value of the CXCL9-CXCR3 axis for Y90-RE-nivolumab combination therapy. Following Y90-RE-nivolumab treatment, R had a higher proportion of hot tumors with increased IFN-γ and cytolytic activity. This change implies successful immune infiltration and activation, which is consistent with the reported mechanisms of action of Y90-RE and the PD-1 inhibitor nivolumab,16 17 47 thereby providing evidence of the conversion of immunologically cold tumors to hot tumors by Y90-RE-nivolumab combination therapy. In accordance with our finding, Chow et al also demonstrated that the CXCL9-CXCR3 interaction enhances intratumoral CD8+ T-cell responses. While CXCR3-deficient mice initially responded poorly to anti-PD-1 treatment, induction of CXCR3 ligands promoted responsiveness.48 The study further implicated CXCL9-CXCR3 as a biomarker of sensitivity to anti-PD-1 treatment. As a primary source of CXCL9 and in response to IFN-γ, TAMs play an important role in antitumoral immunity by recruiting and positioning various activated immune cells via CXCR3 (Cytotoxic T cells (CTLs), NK cells, macrophages, etc) from the surrounding near antigen-presenting cells.49 Furthermore, CXCL9 may also provide additional signals to CTL to facilitate activation and increased IFN-γ secretion. This continuous cycle of paracrine signaling enhances immune cell activation, infiltration, and IFN-γ and CXCL9 secretion, thereby promoting antitumoral responses.50 Taken together, these findings suggest the value of the CXCL9-CXCR3 interaction and macrophages for predicting responses to Y90-RE-nivolumab combination therapy.

In our analysis of the on-treatment biopsy samples at both the RNA and protein levels, we did not identify any immune lineages associated with treatment response, possibly driven by the small sample size as well as extremely hot immune landscape in the tissue. However, the blood-based immune markers indicated the potential of an on-treatment liquid biopsy approach for distinguishing R from NR, with the timing of blood collection matched to that of the tissue biopsy. For example, CXCL11 is often significantly upregulated, which in turn promotes the proliferation and drug resistance of HCC tumor-initiating cells.51 52 Good prognosis in patients with HCC is associated with pro-inflammatory T-helper 1 cytokines, such as IL-12p70 and IL-1β, whereas poor prognosis is associated with T-helper 2 cytokines.53 Even though our data did not confirm all the biomarkers from previous studies, for example, baseline levels of IP-10 (CXCL10) was associated with a better survival in patients with HCC treated with Y90,54 it was found to be non-significant in our study, the existence of many significant biomarkers in the blood highlight the promise of blood-based immunomonitoring as a less invasive method for evaluating responsiveness to the combination treatment. Furthermore, given the challenges associated with obtaining on-treatment tissue biopsies in the real-world setting, this might be a good news for oncologists as well as patients.

By combining analysis of biomarkers across multiple layers, we observed distinct differences in genotypic, phenotypic as well as microenvironmental levels between responders and non-responders (figure 6C). Even though the sample size is not large, the significance of these biomarkers is relatively robust even when we control for confounding variables such as extrahepatic spread (online supplemental figure 15). Moreover, many of these biomarkers had prognostic abilities (online supplemental figures 16 and 17) and tend to act complementarily with each other (eg, TMB and NCOR1 mutations are positive predictors, but chromosome 16 deletions are negative predictors). Among all the biomarkers, the chromosome 16 deletion and the Kaya_P2 subtype had a strong negative correlation with the treatment response. Previous comparison between Asian and European HCCs has revealed a suite of genomic changes associated with this transcriptomic subtype, including higher AFP, AXIN1 mutations as well as higher frequencies of MDSC (Myeloid-Derived Suppressor Cells) Derived suppressor cells) .24 As this subgroup of patients are immunologically cold and might be refractory to current immunotherapy-based modalities, it will be important to look for alternative therapeutic options for these patients. After excluding patients with a high Kaya_P2 signature in our cohort, we observed a response rate of 54% (7 responders out of 13 patients, figure 6C), suggesting possible synergy between Y90-RE and nivolumab in a subset of patients. Furthermore, the response rate can be even higher if we excluded patients with extrahepatic volume (figure 6C). Given the strong ethnic differentiation in HCC24 and its potential link to ICB (Immune Checkpoint Blockage) response,55 the observations presented here suggest great opportunities for tailoring the combination treatment to patients from different ethnic backgrounds.

There are a number of limitations to the current study worth discussing here. First, even though many biomarkers were discovered in this work, the sample size is rather modest and patient variety is limited (eg, patients were recruited from a single center from Singapore with a preponderance of viral positive cancers). As one of our colleagues from Spain is conducting a similar trial, but with a slightly different patient composition (eg, viral negative HCCs),15 a future integrative study combining both cohorts will be able to confirm many of these findings and understand the relative importance of different predictors. Second, even though one of the strengths of this study is paired pretreatment and on-treatment samples and we did observe a convergent response in this cohort to patients treated with Y90 only (eg, elevated CD8+ T cells as well as CCL5 and CXCL16 expression, (online supplemental figure 18).16 However, as the on-treatment samples were collected after the combination therapy, without proper control (ie, on treatment tissue samples after Y90), it is still quite difficult to disentangle the treatment effects and understand whether the observed changes are due to radiation or immunotherapy. Third, even though the significant biomarkers provide several interesting hypotheses about the favorable treatment response, completely understanding molecular mechanisms behind the favorable response will require further studies. For example, radioembolization can lead to activation of immune system (eg, higher expression of interferon pathways) and ‘leaking’ of neoantigens (eg, rapid death of cancer cells), tumors with higher immune activation or higher TMB are more likely to have a better treatment response. In addition, NCOR1 mutations seem to point to a distinctive subset of patients whose response mechanisms might involve complex interplay between the TME and the combination therapy.21

Despite the modest size of the current cohort, we discovered many genomic changes that correlated strongly with treatment response. It is noteworthy that the predictive biomarkers identified for sequential Y90-RE and nivolumab in this study differ from those derived from earlier immunotherapy.20 31–33 For example, high expression of PD-L1, a T-effector signature, and intratumoral CD8+ T-cell density were associated with better outcomes following combination therapy with atezolizumab and bevacizumab.19 In addition, responders to nivolumab and cabozantinib combination therapy were defined by enrichment in effector T cells, tertiary lymphoid structures, CD138+ plasma cells, and a distinct spatial arrangement of B cells.56 However, in our cohort, only the CXCL9-CXCR3 axis and CD68+ macrophages demonstrated the potential for distinguishing responders from non-responders, while neither PDL1-PD-1 axis nor CD8+ T cells offered any predictive value. Furthermore, TMB has not been associated with treatment response in previous immunotherapy-based trials in HCC,19 but showed predictive ability in this study. In addition, the predictive baseline gene signatures known for HCC and other cancer types are not significant in the current cohort (figure 3H). These findings suggest that while common tumor-agnostic biomarkers have been identified, indication-specific and treatment-specific biomarkers are expected. Hence, given the increasing variety of therapeutic options available, a personalized biomarker-directed treatment approach should be our aspiration.


Patient recruitment and sequencing

The clinical trial (CA 209–678) was a single-arm, single-center, two-stage phase 2 trial designed to investigate the activity and safety of Y90-radioembolization followed by nivolumab in patients with advanced HCC. Participants provided written informed consent.17 Thirty-six patients with Child-Pugh A cirrhosis and advanced HCC not suitable for curative surgery were recruited from the National Cancer Centre Singapore/Singapore General Hospital, Singapore. The full inclusion/exclusion criteria and response determination based on the RECIST method can be found in the original study report.17 For each patient, pretreatment and on-treatment biopsies were collected from the same Y90-exposed liver lesion before Y90-RE and after one dose of intravenous nivolumab (figure 1A). AllPrep DNA/RNA Mini Kits (Qiagen) were used to extract DNA and mRNA from tissue biopsies. The NimbleGen SeqCap EZ Human Exome Library (V.3.0) was used for exome capture. Sequencing libraries were prepared using the sequencing platform at the Genome Institute of Singapore and sequenced using the Illumina HiSeq 4000 system.

Somatic variant calling

Short sequencing reads were mapped to the human reference genome (hg19) using bwa.57 Somatic single nucleotide variants were called using Mutect V.1.258 and indels were called using Strelka.59 Somatic variants were subsequently annotated using Oncotator.60 TMB was calculated by dividing the total number of non-synonymous coding mutations by the size of the coding region (ie, 30 Mb).

Neoantigen prediction

Class I HLA (Human Leukocyte Antigens) alleles of patients were identified using Polysolver from patients with genomic information.61 Neoantigens were predicted from somatic mutations based on the inferred HLA types using NetMHCpan (V.4.0) with peptide sizes in the range of 8–11. Neoantigens with a binding affinity of less than 500 nm were kept for further analysis.62 Neoantigen burden was calculated as the total number of unique single nucleotide mutations with predicted neoantigens.

Copy number alterations and genome instability metrics

Copy number alterations (CNA) were identified using Sequenza.63 The GII was calculated as the fraction of the genome having a copy number profile different from the median ploidy across the genome. SCNA scores were calculated for chromosomal arm and broad scale CNV output with the GISTIC algorithm64 using a method similar to that described by Yuan et al.65

CCF and mutational classes under clonal evolution

CCF was calculated using PyClone-VI.66 To delineate the somatic variants into persistent, expanding and contracting mutations, we subtracted the pretreatment CCF of each mutation from its corresponding on-treatment CCF value. If the difference was <10% of the baseline value, the mutation was categorized as a persistent mutation. If the CCF increased by ≥10%, the mutation was categorized as an expanding mutation. If the CCF was decreased by ≥10%, the mutation was categorized as a contracting mutation.

RNA-seq analysis

Sequencing reads were mapped to the hg19 reference genome using STAR67 and quantified using RSEM.68 DEGs were identified using DESeq2 (log2 fold-change >1 or <−1 and adjusted p value<0.05). Pathway enrichment analysis was conducted using the ‘fgsea’ R package69 with HALLMARK gene sets from the MSigDB database.70 Signatures for previously reported HCC subtypes24 26 27 were downloaded from the MSigDB database. Gene signatures for subtypes in Kaya et al24 were obtained as top 100 upregulated genes when comparing each subtype (eg, P1, P2, M1, M2) with the rest of the subtypes in Asian patients of TCGA-LIHC cohort (online supplemental table 2). Previously reported immunotherapy and radiotherapy gene signatures were collected from various studies, including signatures reported by Sangro et al (n=7),32 using the TIDE tool online (n=11),30 two publications on radiosensitivity scores (n=2),33 34 gene expression profiling score71 and a publication reporting a decrease in a nine-gene WNT score following immunotherapy.31 All gene signature scores were calculated by gene set variation analysis (GSVA) using the ‘gsva’ R package.72 Bulk RNA-seq deconvolution to 15 immune cell types was conducted using a method described in Danaher et al.36

Flow cytometry

Peripheral blood mononuclear cells (PBMCs) were isolated by the Ficoll-Paque density gradient centrifugation method and stained as previously described.73 Briefly, samples were incubated with Zombie NIR Fixable Viability dye (BioLegend, San Diego, California, USA) for 10 min at 4°C in the dark for live/dead cell discrimination. Fc receptors were blocked with Human TruStain FcX (BioLegend) for 10 min at room temperature. Cell surfaces were labeled with antibodies targeting markers of interest (online supplemental table 3) for 30 min at 4°C. Samples were fixed and permeabilized using BD Cytofix/Cytoperm (BD Biosciences, Franklin Lakes, New Jersey, USA) before staining with antibodies targeting intracellular markers (online supplemental table 3). Data were acquired using the Cytek Aurora spectral flow cytometer (Cytek Biosciences, Fremont, California, USA) and analyzed using FlowJo V.10 software (FlowJo LLC, Ashland, Oregon, USA).

Multiplex immunohistochemistry

The mIHC was performed using an Opal Multiplex fIHC kit (Akoya Biosciences, Marlborough, Massachusetts, USA) on a Leica Bond Max autostainer (Leica Biosystems, Melbourne, Australia) as previously described.73–76 In brief, deparaffinized/rehydrated formalin-fixed, paraffin-embedded tissue sections were subjected to heat-induced epitope retrieval, peroxidase blocking, and incubation with primary antibodies against CD8, CD38, CD45, CD68, CXCR3, CXCL9, LAG-3, PD-1, PD-L1, or STAT1 (online supplemental table 4), followed by the addition of polymeric horseradish peroxidase-conjugated secondary antibodies (Leica Biosystems, Newcastle-upon-Tyne, UK) and Opal tyramide signal amplification (TSA) reagents (Akoya Biosciences). Following TSA deposition, the slides were then subjected to heat-induced stripping and the labeling processes were repeated until six markers were labeled prior to counterstaining with spectral DAPI (Akoya Biosciences). Slides were visualized using a Vectra 3 pathology imaging system (Akoya Biosciences). A pathologist (JY) assessed the slides visually and selected multiple regions of interest (ROI) with high-quality staining and viable tumor cells. These ROIs were scanned at 20× magnification and then analyzed and scored by the pathologist using inForm software (V.2.4.2; Akoya Biosciences) and HALO (Indica Lab, Albuquerque, New Mexico, USA).

Multiplex cytokine analysis

Plasma was collected during the PBMC isolation process as previously described.77 Briefly, whole blood was centrifuged for 15 min at 400×g and the upper (plasma) layer was harvested. Following centrifugation at 1,000×g for 10 min, the supernatant was collected and stored at −80°C prior to further analysis. A total of 65 analytes from the Immune Monitoring ProcartPlex Panel (Thermo Fisher Scientific, Waltham, Massachusetts, USA) were analyzed according to the manufacturer’s protocol (online supplemental table 5). After plates were washed using a Tecan Hydrospeed Washer (Tecan Group AG, Männedorf, Switzerland), data were acquired with a Flexmap 3D system (Luminex Corporation, Austin, Texas, USA) and analyzed using Bio-Plex manager V.6.2 software (Bio-Rad Laboratories, Hercules, California, USA) with a 5-parameter curve-fitting algorithm applied for standard curve calculations.

Statistical analysis

For patients without paired biopsies, a partially overlapping t-test was applied when comparing pretreatment and on-treatment data for gene signatures and genomics features. For unpaired comparison, p values from Wilcoxon rank-sum test are shown on the boxplots. Fisher’s exact test was implemented for two categorical features. When there are multiple samples from one patient at a single treatment point, average value was used for comparisons.

Data availability statement

Data are available in a public, open access repository. Raw whole exome sequencing (WES) and RNA sequencing fastq files are deposited in The European Genome-phenome Archive (EGA) under the study with accession code ‘EGAS00001006834’.

Ethics statements

Patient consent for publication

Ethics approval

This study involves human participants and was approved by SingHealth Centralised Institutional Review Board IRB Ref. No. 2016/2613. Participants gave informed consent to participate in the study before taking part.


We thank Anders Jacobsen Skanderup and Qianfei Wang for their constructive feedback and discussions. We acknowledge the support from Bristol Myers Squibb and Sirtex.


Supplementary materials

  • Supplementary Data

    This web only file has been produced by the BMJ Publishing Group from an electronic file supplied by the author(s) and has not been edited for content.


  • NAK, DT and XL contributed equally.

  • Contributors Conceptualization: NAK, DT, JY and WZ. Data collection and analysis: NAK, DT, XL. JQL, JY and WZ. Interpretation of results: NAK, DT, XL, JQL, MCL, DG, CZJP, FW, CRJ, JCTL, ZWN, JY, CCLC, JJXL, KSHL, AG, FY, MYL, TWHS, HCT, AH, YSC, TK-HL, WLT, SPC, JY and WZ. Writing—Original Draft: NAK, DT, XL, JY and WZ. Writing—Review and Editing: NAK, DT, XL, DG, CCLC, CZJP, MYL, AH, YSC, WLT, SPC, JY and WZ. Supervision: DT, JY and WZ. WZ is the overall guarantor for this article.

  • Funding This study is funded by the National Medical Research Council grant NMRC/CIRG/1454/2016. WZ is supported in part by the National Natural Science Foundation of China (32293190/32293192, 31970566), the Strategic Priority Research Program of the Chinese Academy of Sciences (XDPB17), National Key R&D Program of China (2018YFC1406902). NAK is supported in part by the Singapore National Medical Research Council grants (TCR/015-NCC/2016). JY is supported in part by Singapore National Medical Research Council grants (MOH-000323-00, OFYIRG19may-0007, NMRC/OFLCG/003/2018) and A*STAR Biomedical Research Council IAF-PP grant (HBMS) (H19/01/a0/024 – SInGapore ImmuNogrAm for Immuno-Oncology (SIGNAL)).

  • Competing interests DT reports receiving research support from Novartis, Sirtex, and Bristol Myers Squibb. DT reports having received honorarium as an advisor for Novartis, Celgene, Sirtex, Merck Sharp & Dohme, Eisai, Ipsen, Bayer, and Bristol Myers Squibb.

  • Provenance and peer review Not commissioned; externally peer reviewed.

  • Supplemental material This content has been supplied by the author(s). It has not been vetted by BMJ Publishing Group Limited (BMJ) and may not have been peer-reviewed. Any opinions or recommendations discussed are solely those of the author(s) and are not endorsed by BMJ. BMJ disclaims all liability and responsibility arising from any reliance placed on the content. Where the content includes any translated material, BMJ does not warrant the accuracy and reliability of the translations (including but not limited to local regulations, clinical guidelines, terminology, drug names and drug dosages), and is not responsible for any error and/or omissions arising from translation and adaptation or otherwise.