Single cell RNA-seq analysis identifies a noncoding RNA mediating resistance to sorafenib treatment in HCC

© The Author(s) 2021. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Main text Sorafenib, a multikinase inhibitor, is an FDA-approved first line drug for the treatment of patients with unresectable hepatocellular carcinoma (HCC) [1]. Unfortunately, these patients achieve minimal therapeutic benefit with an improvement in overall survival of only 3 months [2]. After an initial response, the majority of HCC patients develop drug resistance and disease progression [3]. Although multiple mechanisms such as crosstalk involving PI3K/Akt and MAPK/ERK pathways or activation of EMT have been reported [4, 5], the major drivers of sorafenib resistance remain obscure. Recently, rapidly evolving single-cell transcriptome analysis techniques have been used to comprehensively assess the genetics of the tumour microenvironment and the mechanisms for intra-tumoral heterogeneity which can impact anti-cancer treatment responses [6, 7]. This technology has the potential for higher transcriptomic resolution and sensitivity than bulk sequencing techniques [6, 7]. Single-cell analysis also enables a deep interrogation of hundreds of cancer drivers and identifies individual cells or genes with specific genetic modifications or expression profiles which could contribute to the development of therapeutic resistance. Here, we conducted a single cell RNA seq analysis in sorafenib resistant HCC cells to systematically investigate drug-resistance signatures to uncover the mechanisms of sorafenib resistance. Results and discussion We successfully established Hep3B and Huh7 sorafenib resistant cells (Hep3B-R and Huh7-R) following longterm exposure of sorafenib to parental Hep3B and Huh7 (Hep3B-P and Huh7-P) cells. The IC50 of Hep3B-R and Huh7-R cells was nearly 2-fold and 2.3-fold higher than that of Hep3B-P and Huh7-P cells (Fig. 1A and Supplementary Fig. 1A). Interestingly, Hep3B-R and Huh7-R cells had greater self-renewal capacity developing significantly larger and more numerous tumour-spheres than the Huh7-P cells (Supplementary Fig. 1B, C). Putative liver cancer stem-like cell populations (LCSCs; defined as EpCAM+CD133+ cells) were found to be significantly increased in both Hep3B-R and Huh7-R cells (Supplementary Fig. 1D). To assess the enrichment of stem-like phenotypes, nanostring analysis was performed to investigate the targeted mRNA expression of 21 stemness related genes (eg., ALDH1A1, ALDH2, ABCC1, CD133, CD44, EpCAM, Oct4, Sox2, β-catenin). As expected, marked to moderate upregulation of majority of these genes was observed in Hep3B-R and Huh7-R cells (Supplementary Fig. 1E). Consistently, we observed that protein levels of EMT or stemness markers like Vimentin, N-cadherin, Notch2, Notch3, AFP and Oct4 were all upregulated, whereas E-cadherin was downregulated in Hep3B-R and Huh7-R cells (Supplementary Fig. 1F). These data indicate that the HCC cells likely undergo EMT and this is associated with the acquisition of stemness features in response to sorafenib treatment. The enrichment of stemness/EMT phenotypes at the population level prompted us to hypothesize that their induction might be a major contributor to sorafenib Open Access


Main text
Sorafenib, a multikinase inhibitor, is an FDA-approved first line drug for the treatment of patients with unresectable hepatocellular carcinoma (HCC) [1]. Unfortunately, these patients achieve minimal therapeutic benefit with an improvement in overall survival of only 3 months [2]. After an initial response, the majority of HCC patients develop drug resistance and disease progression [3]. Although multiple mechanisms such as crosstalk involving PI3K/Akt and MAPK/ERK pathways or activation of EMT have been reported [4,5], the major drivers of sorafenib resistance remain obscure. Recently, rapidly evolving single-cell transcriptome analysis techniques have been used to comprehensively assess the genetics of the tumour microenvironment and the mechanisms for intra-tumoral heterogeneity which can impact anti-cancer treatment responses [6,7]. This technology has the potential for higher transcriptomic resolution and sensitivity than bulk sequencing techniques [6,7]. Single-cell analysis also enables a deep interrogation of hundreds of cancer drivers and identifies individual cells or genes with specific genetic modifications or expression profiles which could contribute to the development of therapeutic resistance. Here, we conducted a single cell RNA seq analysis in sorafenib resistant HCC cells to systematically investigate drug-resistance signatures to uncover the mechanisms of sorafenib resistance.

Results and discussion
We successfully established Hep3B and Huh7 sorafenib resistant cells (Hep3B-R and Huh7-R) following longterm exposure of sorafenib to parental Hep3B and Huh7 (Hep3B-P and Huh7-P) cells. The IC50 of Hep3B-R and Huh7-R cells was nearly 2-fold and 2.3-fold higher than that of Hep3B-P and Huh7-P cells ( Fig. 1A and Supplementary Fig. 1A). Interestingly, Hep3B-R and Huh7-R cells had greater self-renewal capacity developing significantly larger and more numerous tumour-spheres than the Huh7-P cells ( Supplementary Fig. 1B, C). Putative liver cancer stem-like cell populations (LCSCs; defined as EpCAM + CD133 + cells) were found to be significantly increased in both Hep3B-R and Huh7-R cells (Supplementary Fig. 1D). To assess the enrichment of stem-like phenotypes, nanostring analysis was performed to investigate the targeted mRNA expression of 21 stemness related genes (eg., ALDH1A1, ALDH2, ABCC1, CD133, CD44, EpCAM, Oct4, Sox2, β-catenin). As expected, marked to moderate upregulation of majority of these genes was observed in Hep3B-R and Huh7-R cells (Supplementary Fig. 1E). Consistently, we observed that protein levels of EMT or stemness markers like Vimentin, N-cadherin, Notch2, Notch3, AFP and Oct4 were all upregulated, whereas E-cadherin was downregulated in Hep3B-R and Huh7-R cells ( Supplementary Fig. 1F). These data indicate that the HCC cells likely undergo EMT and this is associated with the acquisition of stemness features in response to sorafenib treatment.
The enrichment of stemness/EMT phenotypes at the population level prompted us to hypothesize that their induction might be a major contributor to sorafenib Open Access resistance. To test this hypothesis, single-cell RNA sequencing was performed in Huh7-P and Huh7-R cells using the 10X Genomic Chromium System. After stringent filtering of low-quality cells, 4776 cells were selected for subsequent single cell sequencing data analysis. The distinct gene expression signatures between Huh7-P and Huh7-R are visualised in two heatmaps ( Fig. 1B, C) and a two-dimensional tSNE plot (Fig. 1F). As shown, Huh7-R cells showed upregulated expression of genes related to stemness and EMT, including multidrug resistant genes (eg., ABCC1, ABBC2, ABBC3), stemness markers (eg., ALDH2, CD13, CD44, CD133, Krt19) and EMT-related genes (eg., LAMB1, LAMC1, S100P, S100A6, S100A11, Vimentin, N-cadherin, ZEB1) (Fig. 1B). It is now understood that Notch signaling is one of the most activated pathways in cancer, playing a central role in cancer progression, EMT promotion and the activation of CSCs [8,9]. Our data also confirmed that the mRNA expression of Notch receptors (Notch2-4), ligands (DLL1, DLL4, Jagged1) and Notch transcriptional genes (Hes1, Hey1) are all upregulated in patient-derived HCCs (n = 12) as compared to matching normal liver, whereas the Notch inhibitor DLL3 was downregulated ( Supplementary   Fig. 2). Huh7-R cells also exhibited enhanced expression of Notch1-3 receptors (Fig. 1D) and various Notch activators including its ligands (Jagged1), Deltex genes (DTX2-4), transcriptional genes (Hes1, Hes4, MAML, ADAM17), GSK3β and γ-secretase complex components (PSEN1, PSEN2) (Fig. 1C). This was accompanied by downregulation of the key Notch repressors DLL3 and RBPJ (Fig. 1C). This indicated that the activation of Notch signaling serves a crucial role in the induction of stemness/EMT traits and the acquisition of sorafenib resistance.
Despite showing upregulation of stemness/EMT genes, Huh7-R cells showed downregulation of the proliferative marker MKI67, and cell-cycle related genes (Fig. 1E) suggesting that sorafenib treatment may enrich a subset of resistant cells with slow cell cycle or those that lie in a dormant state. This is consistent with the phenomenon of quiescent stem-like cells being maintained in a dormant or slow-dividing state, but subsequently expanding and dividing following cancer treatments [10]. To further confirm the presence of these cells, we identified a small fraction of Huh7-P cells (as indicated by black arrows) which were clustered together with the population of Arrows indicate a few Huh7-P cells (candidate stem cells) that are clustered together with Huh7-R cells. G Dot-plot showing the expression levels of stemness/EMT genes in Huh7-P candidate stem cells and Huh7-R cells Huh7-R cells (Fig. 1F). We next examined the expression of stemness/EMT genes in these individual cells. Surprisingly, these rare Huh7-P cells (named as candidate stem cells) showed gene expression signatures with great similarity to that of Huh7-R cells (Fig. 1G) suggesting that drug-resistant cells with stemness/EMT traits may have pre-existed in the naïve cells and become prevalent following long-term sorafenib exposure.
Among the genes related to sorafenib resistance, we focused on ZNFX1 antisense RNA 1 (ZFAS1) which had the highest upregulation in Huh7-R cells. ZFAS1 as a novel regulator IncRNA was found to be overexpressed in HCC (n = 10) (Supplementary Fig. 3A) and is reported to be involved in tumour growth and metastasis in multiple cancers [11,12]. We therefore determined the correlation of ZFAS1 expression and clinicopathological features in HCC by analysing publicly available datasets. In this analysis, we found that ZFAS1 expression was significantly upregulated in HCC and recurrent tumours compared with that of matching normal liver by using the GEPIA dataset ( Supplementary Fig. 3B) and the University of California Santa Cruz (UCSC) Xena project (http:// xena. ucsc. edu) (Supplementary Fig. 3C). Moreover, ZFAS1 was overexpressed at higher stages (III/IV) of HCC and in poorly (G3)/undifferentiated HCC (G4) compared with those at early stages (I/II) (Supplementary Fig. 3D) and in highly (G1)/moderately (G2) differentiated HCC (Supplementary Fig. 3E). By analysing the GEPIA dataset, higher ZFAS1 expression levels were associated with shorter overall survival (OS, Supplementary Fig. 3F) and disease-free survival (DFS, Supplementary Fig. 3G) suggesting that ZFAS1 may serve as a novel prognostic marker of HCC. To clarify the role of ZFAS1 in the induction of stemness and EMT phenotypes, we used the cBioportal dataset to investigate the correlation of ZFAS1 and stem-ness/EMT genes. As expected, ZFAS1 expression was positively correlated with the expression of multiple stemness genes (eg., EpCAM, CD24, CD90, CD133, DLK1, Krt18/19) and EMT markers (eg., S100A4/A6, Twist, Vimentin) ( Supplementary Fig. 4).

Conclusions
In summary, this study has demonstrated the precise role of stemness and EMT in driving sorafenib resistance at the single cell level. A small fraction of quiescent stemlike cells with intrinsic sorafenib resistance pre-exists in the HCC cell population and becomes enriched following long-term sorafenib exposure. This data provides new insights indicating that dormant stem-like cells are possibly the main culprits facilitating the development of sorafenib resistance. Therefore, killing quiescent cancer stem-like cells might be a strategy to overcome sorafenib resistance. We also identified a novel noncoding lncRNA, ZFAS1 which is closely associated with HCC proliferation and the maintenance of stemness characteristics. Silencing ZFAS1 is thus another strategy to overcome sorafenib resistance.