Immunosuppression and phenotypic plasticity in an atlas of human hepatocholangiocarcinoma

Introduction

Primary hepatic carcinoma is one of the most common cancers worldwide (1), and hepatocholangiocarcinoma (H-ChC) is a specific form of hepatic carcinoma exhibiting the pathological features of both hepatocellular carcinoma (HCC) and intrahepatic cholangiocarcinoma (iCCA) (2). The true incidence of H-ChC may be underestimated because it may be misdiagnosed as HCC or iCCA, and most patients are inoperable (3). Moreover, the specific risk factors associated with its development remain unknown. Some studies have shown that the expression of nestin, a marker of dual-potent progenitor oval cells, can be used as a diagnostic biomarker of H-ChC (4,5). However, nestin is expressed in only a subset of HCC and iCCA cases, and its value for H-ChC remains to be evaluated. Therefore, the molecular characterization of H-ChC is urgently needed to facilitate accurate diagnosis and accelerate the development of professional therapy.

The clonal origin of H-ChC remains a matter of debate. H-ChC tumours exhibit a high degree of intratumoural heterogeneity, and features associated with different subtypes often coexist within the same tumour; further studies are needed to better identify the different expansions associated with different phenotypes. Although genetic alterations may drive histological differentiation, a role of the microenvironment (inflammatory and immune cells, fibroblasts, and blood vessels) cannot be ruled out. Single-cell RNA sequencing (scRNA-seq) can provide insight into these issues by revealing the intratumour heterogeneity and tumour immune microenvironment (TIME) of H-ChC.

In this study, we performed scRNA-seq analysis using tissue specimens which got from the Department of Liver Surgery of Peking Union Medical College Hospital (PUMCH), obtaining H-ChC patients and tissue samples of HCC and iCCA patients used as controls. We constructed a comprehensive atlas for H-ChC, in which the mixed lineage tumour cells of three tumours were detected for the first time. Our data indicate two subtypes of H-ChC: one has a metabolism-related phenotype (CHP), and the other has a proliferation-related phenotype (CIP), suggesting that H-ChC could be identified through CIP and CHP at the whole-transcriptome level. In addition, a specific tumour cell subpopulation enriched in CHP has characteristics of activation status [named non-quiescent cell cycling malignant epithelial cells (NCM-Epi)], which is more strongly related to disorders of CD8⁺ Tex cells in H-ChC. We also found that GZMK⁺ CD8⁺ Tex cells may play a key role in clinical immunotherapy, and transcriptome heterogeneity might be largely driven by NCM-Epi and GZMK⁺ CD8⁺ Tex cell dysfunction. Our analysis is critical for understanding the carcinogenesis of H-ChC and presents a thorough characterization of its molecular features.

Methods

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Haitao Zhao (pumchzht@aliyun.com).

Material availability

This study did not generate unique reagents.

Experimental model and subject details

Human specimens

In the Department of Liver Surgery at PUMCH, 13 cancer patients were enrolled and pathologically diagnosed with H-ChC, HCC, or iCCA. The patients included 11 males and 2 females. Fresh tumour samples were surgically dissected from the aforementioned patients. A median age of 57 years was observed among the participants, and patient age ranged from 37 to 78 years. According to the American Joint Committee on Cancer (AJCC) and Barcelona Clinic Liver Cancer (BCLC) recommendations, individuals with HCC and iCCA were divided into distinct clinical phases. Patients with unresectable advanced H-ChC confirmed by pathology were included in the study, while patients with HCC or cholangiocarcinoma were excluded. Patients who had received chemotherapy or radiation before tumour removal were also excluded from the study. The 10 patients with advanced H-ChC in the cohort study had received immunotherapy-based systemic therapy during treatment and had more than 1 measurable lesion [according to Response Evaluation Criteria In Solid Tumors, version 1.1 (RECIST v1.1) criteria]. Table S1 summarizes the clinical features and follow-up information. The studies involving human participants were reviewed and approved by the institutional review board and ethics committee at Peking Union Medical College Hospital (No. JS-1391). The patients/participants provided their written informed consent to participate in this study. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Method details

Single-cell collection

Tissue samples were obtained by washing in ice-cold RPMI 1640 and dissociating using a Multitissue Dissociation Kit 2 following the package recommendations. The viscosity of the homogenate was taken into consideration during DNase treatment. After removing erythrocytes from the samples (Miltenyi 130-094-183, Shanghai, China), the cells were counted, and viability was determined using a fluorescence cell analyser (Countstar^® Rigel S2, Shanghai, China) with acridine orange (AO)/propidium iodide (PI) reagent. Removal of debris and dead cells was then performed as needed (Miltenyi 130-109-398/130-090-101). Finally, after a second wash in RPMI 1640, fresh cells were resuspended in 1× PBS containing 0.04% bovine serum albumin at 1×10⁶ cells per mL.

Generation of patient scRNA-seq data

scRNA-seq libraries were prepared with a SeekOne^® MM Single Cell 3' Library Preparation Kit (SeekGene, Beijing, China). Briefly, an appropriate number of cells were loaded into the flow channel of the SeekOne^® MM chip with 170,000 microwells and allowed to settle in the microwells by gravity. After removing unsettled cells, sufficient cell-barcoded magnetic beads (CBBs) were permitted to stay after being pipetted into the flow path in the microwells with the aid of a magnetic field. Excess CBBs were rinsed out, and the cells in the MM chip were lysed to release RNA, which was captured by the CBBs in the same microwell. Then, all CBBs were collected, and reverse transcription was performed at 37 °C for 30 min to label cDNA with cell barcodes on the beads. Exonuclease I treatment was performed to remove unused primers on the CBBs. Subsequently, barcoded cDNA on the CBBs was hybridized with a random primer carrying the SeqPrimer sequence on the 5' end, which was extended to form the second strand DNA with the cell barcode on the 3' end. The resulting second-strand DNAs were denatured off the CBBs, purified, and amplified by PCR. The amplified cDNA product was then cleaned to remove unwanted fragments, a full-length sequencing adapter was added, and the sample was indexed by index PCR. We cleaned the indexed sequencing libraries using SPRI beads and calculated them using quantitative PCR (KAPA Biosystems KK4824) and an Illumina NovaSeq 6000 with a PE150 read length.

scRNA-seq data processing

SeekOne Tools were used with the scRNA-seq data to match reads and produce gene-cell-specific unique molecular identifiers (UMIs). We utilized the GRCH38 genome as a guide. The Seurat methodology was then applied in subsequent analyses (6). Using the Harmony function, various samples were combined. Cells were subjected to additional quality control; the high-quality cell numbers in each group are listed in Table S1. Using the NormalizeData function, a logarithmic data matrix was generated, and the FindVariableFeatures function was applied to discover genes with high variability. The dimension numbers were reduced by T-distribution stochastic neighbour embedding (t-SNE) implemented with the DimPlot function with default parameters for better visibility (7).

Annotation of cell types

The expression of recognized markers was used to annotate cell types: ALB, TF, TTR, HNF4A, CYP2A6, EPCAM, KRT19 and CD24 for detecting hepatocholagiocyte; ALB, TF, TTR, HNF4A, CYP2A6 for hepatocyte; EPCAM, KRT19 and CD24 for cholangiocyte; PDGFRB, ACTA2, COL1A1, COL1A2, COL3A1, DES and DCN for fibroblasts; PECAM1, CDH5, ICAM2, KDR and ERG for endothelial cells (ECs); CD3D, CD3E and CD3G for T cells; CD79A, CD79B, CD19 and MS4A1 for B cells; CD14, CD163, CD68 and CSF1R for myeloid cells; and KIT, TPSAB1 and TPSB2 for mast cells. The method for identifying the three epithelial cell types is based on the study by Li et al. (8). Cancer cells were identified as aneuploid cells using the R package CopyKAT (9).

Gene programs in malignant epithelial cells

The utilization of nonnegative matrix factorization (NMF) on tumour cells was in contradictory to previous published articles (10), which used NMF independently for every patient, and then used hierarchical clustering to identify consensus expression programs across patients (11). We believed that if we could find that heterogeneity in patient was smaller than heterogeneity in cancer, then it would be more necessary to perform NMF among cancers rather than integrating the results of NMF for each patient. Therefore, we employed two measures, mutual information (MI) and heterogeneity score (12,13), to validate that heterogeneity of patient was smaller than heterogeneity in cancer types. The specific procedure involved: (I) randomly sampling 75% of the cells for 200 rounds to identify characteristic genes; (II) calculate cosine distances (III) NMF was then performed to extract expression programs; (IV) enrichment analysis (14). Then we identified two subtypes CHP and CIP.

Detecting the heterogeneity of malignant cells

We utilized the diversity score to quantify heterogeneity between CHP and CIP (15).

Differences in the proportions of cell subpopulations

We used the scDC algorithm (16) to execute differential cell type composition analysis, and a bias-free and accelerated bootstrap calculation of cell-sort percentage comparison was applied.

Construction of single-cell trajectories

Malignant cell transitions from one state to another were inferred by performing partition-based graph abstraction (Slingshot) analysis (17).

Analysis of cell-cell interactions

The correlation between immune cells and stromal cells was predicted with iTALK (version 2.0) (18) and scCellChat (19).

Identification of signature genes

The marker genes for every group were identified using the Seurat package’s FindAllMarkers function.

Survival analysis

We used the R packages survival and survminer for survival analysis. The R functions survfit() and survdiff() were utilized to produce Kaplan-Meier survival curves, calculate P values, and conduct log-rank testing. The tumour response was evaluated by two independent professional radiologists who were blinded to the therapeutic outcomes and clinicopathological features at PUMCH according to RECIST V1.1. Objective response rate (ORR), disease control rate (DCR), and progression-free survival (PFS) are acknowledged as useful clinical metrics of therapeutic efficacy. ORR is defined as the proportion of patients who achieved objective response, including complete remission (CR) or partial remission (PR); DCR is the proportion of patients who acquired objective response or stable disease (SD); and PFS is the time from the initial dose to the first radiologically confirmed progressive disease (PD) or death. All of the outcomes are determined by radiology.

Quantification and statistical analysis

Statistical analysis

Statistical analyses were conducted as described in the figure legends. The R package pwr was used to perform power analyses.

Results

A comprehensive atlas of different cell components of H-ChC, HCC, and iCCA

The coexistence of HCC and iCCA components in H-ChC is an interesting phenomenon, and scRNA-seq was performed to investigate the relationship between these two components. Samples were obtained from 13 patients in the Department of Liver Surgery of PUMCH (Figure 1A and Table S1), including 8 samples from 5 H-ChC patients, 4 samples from 4 HCC patients, and 4 samples from 4 iCCA patients. We first analysed the 15 samples from the 12 untreated tumour patients and obtained 35 cell subpopulations of 91,112 high-quality single cells, which included 15,708 epithelial cells, 63,248 immune cells, and 12,156 stromal cells (Figure 1B and Figure S1).

Figure 1 A comprehensive atlas of different cell components of H-ChC, HCC, and iCCA. (A) Study details. The clinical characteristics of the sample are shown on the left. Samples from different tumour types were collected. (B) t-SNE plot of all cells coloured by cell type. (C) Histopathology images of samples from H-ChC, HCC, and iCCA patients. Scale bars, 50 µm. (D) The columns show the expression levels of CK19 and Arg-1 in H-ChC, HCC, and iCCA. BCLC, Barcelona Clinic Liver Cancer; HBV, hepatitis B virus; HCV, hepatitis C virus; ECOG, Eastern Cooperative Oncology Group; H-ChC, hepatocholangiocarcinoma; HCC, hepatocellular carcinoma; iCCA, intrahepatic cholangiocarcinoma; IITH, inter- and intra-tumor heterogeneity; CNV, copy number variant; TIME, tumour immune microenvironment; Tex, exhausted T cell; EC, endothelial cell; DC, dendritic cell; t-SNE, T-distribution stochastic neighbour embedding.

The cells were annotated according to known cell-specific marker genes (20-23) of epithelial cells, T cells, B cells, fibroblasts, ECs, and myeloid cells (Figures S1B,S2A). We observed that compared with nonepithelial cells, epithelial cells from different patients formed farraginous clusters and exhibited more apparent inter- and intra-tumor heterogeneity (IITH) (Figure S1A). Inter-tumour heterogeneity refers to the heterogeneity between different tumours, while intra-tumor heterogeneity refers to the heterogeneity among different samples/cells within the same tumour. The same type of cancer may present differently in different patients, and there may also be variations among different tumour sites within the same patient. Subsequently, we observed differences in the proportions of epithelial to nonepithelial cells among the three tumours (Figure S2B-S2F), highlighting the complex composition of the TIME in H-ChC, HCC, and iCCA. Three tumours were verified by histopathological analysis (Figure 1C,1D). All evidence supported the heterogeneous landscape of primary liver cancer (PLC) biopsies observed in the pathology images.

scRNA-seq detected farraginous-lineage tumour cells in H-ChC, HCC and iCCA

To characterize the underlying molecular features of tumour cells exhaustively, we constructed an epithelial cell atlas (Figure 2A-2C) that contained 5,234 hepatocholangiocytes (simultaneously expressing hepatocyte and cholangiocyte markers), 6,268 hepatocytes and 3,527 cholangiocytes (Figure S3A,S3B, Table S2) (15,24-26). At the tumour cell level, to distinguish different pathological features of H-ChC, single-cell copy number variants (CNVs) were acquired through CopyKAT (9), and malignant epithelial cells were then obtained (Figure S3C) (4,5). We found that the proportion of tumour cells we redefined expressing well-established marker genes of hepatocytes and cholangiocytes ranged from 15% to 93.9% depending on the gene expression pattern (Figure 2B-2D, Figure S3D-S3F). Hepatocholangiocytes were observed in all the patients from whom we isolated epithelial cells (Figure 2D), hepatocytes were enriched in HCC, while cholangiocytes were enriched in iCCA. In addition, these results were also confirmed by the previously published scRNA-seq dataset (27), in which we used the same method to calculate the expression of the cancer-specific marker genes. A total of 69.47% (2,301 out of 3,312) of the tumour cells in H-ChC were considered hepatocholangiocytes (Figure S3G-S3I). To investigate this further, the expression of marker genes of hepatocytes and cholangiocytes were compared among the three cell types (hepatocytes, cholangiocytes, and hepatocholangiocytes). We have demonstrated through RNA velocity that hepatocholangiocytes are a mixed cell type between hepatocytes and cholangiocytes (Figure 2E). Using previously established methods, we constructed a Markov graph from the dataset. The expression pattern of the corresponding gene can be readily represented by a 3-coordinate ternary graph, as implemented in the ggtern package (Figure 2F) (12). To the best of our knowledge, this is the first work to detect farraginous-lineage tumour cells at the single-cell level in liver cancer. All the results above showed that hepatocholangiocytes are not unique in H-ChC but are prevalent in most PLC patients, although not in the same proportion.

Figure 2 The farraginous-lineage tumour cells in H-ChC, HCC and iCCA. (A) t-SNE plot of epithelial cells coloured according to cell type, hepatocholangiocytes, hepatocytes, and cholangiocytes. (B) The proportions of epithelial cell types in different patients. Each bar is coloured according to cell type. (C) Bubble heatmap showing the value of hepatocyte score and cholangiocyte score in distinct epithelial cell subtypes. (D) Boxplots showing the hepatocyte score and cholangiocyte score alterations of tumour cell clusters in different specimens. *, P<0.05; **, P<0.01; ***, P<0.001. (E) Hepatocholangiocytes tend to be a mixed cell type between hepatocytes and cholangiocytes, as shown by CellRank. (F) The expression pattern of the corresponding gene in hepatocholangiocytes, hepatocytes and cholangiocytes. H-ChC, hepatocholangiocarcinoma; HCC, hepatocellular carcinoma; iCCA, intrahepatic cholangiocarcinoma; t-SNE, T-distribution stochastic neighbour embedding.

H-ChC has two phenotypes: CHP and CIP

Subsequently, to better identify the effect of the malignant epithelial component on H-ChC tumour progression, the status of the 9,347 malignant epithelial cells (hepatocholangiocytes: 4,372; hepatocytes: 2,271; cholangiocytes: 2,424) was characterized at single-cell resolution (Figure S4A). It has been suggested that tumour cells cannot be used directly for stratification (27), and we observed that the tumour cells of HCC, iCCA, and H-ChC were not divided by cancer subtype due to the existence of hepatocholangiocytes (Figures S3D,S4A-S4K), although the heterogeneity between patients is significantly lower than that between cancer types (Figure S5A-S5C).

Because the single-gene approach is often unreliable for scRNA-seq due to common gene deletions (22,23,28,29), the cellular composition of different components and gene programs of tumour cells were assessed by NMF. The NMF results showed (I) an H-ChC component that is functionally similar to HCC with high expression of iCCA epithelial markers (i.e., an HCC component with the iCCA phenotype, CHP); (II) an H-ChC component that is functionally similar to iCCA with high levels of HCC epithelial markers (i.e., an iCCA component with the HCC phenotype, CIP) (Figure 3A,3B, Figure S4B,S4C); and (III) CHP and CIP were both enriched in H-ChC, while HCC just comprised CHP, iCCA contained CIP merely. The relationship between the hepatocyte/cholangiocyte scores and the CHP/CIP phenotype, reflecting the distribution of the three cell types in CHP and CIP based on the two scoring measures (Figure 3B).

Figure 3 Two phenotypes of H-ChC: CHP and CIP. (A) Two gene programs (CHP and CIP) classified by NMF. (B) The relationship between the hepatocyte/cholangiocyte scores and the CHP/CIP phenotype, reflecting the distribution of the three cell types in CHP and CIP. (C) The enriched gene pathways of CHP and CIP. (D) Intratumoural heterogeneity within different malignant cell types as measured by diversity score for CHP and CIP. (E) Differential expression analysis of two phenotypes in different clusters. (F) Boxplots showing the patterns of STMN1, EMT and cell cycle score in different cell clusters of malignant epithelial cells. (G) Heatmap showing the 5 significantly activated transcription factors in NCM-Epi. (H) TF target network created from NCM-Epi and other malignant cells, organized into super-regulons for NCM-Epi. H-ChC, hepatocholangiocarcinoma; HCC, hepatocellular carcinoma; iCCA, intrahepatic cholangiocarcinoma; NMF, nonnegative matrix factorization; EMT, epithelial mesenchymal transition; NCM-Epi, non-quiescent cell cycling malignant epithelial cells.

In this context, the biphenotype refers to the functional similarity while having high expression of iCCA/HCC marker genes in the same tumour cell population, rather than being an admixture of two or more tumour cell subpopulations. The findings support the point of view that H-ChC can be defined as a separate tumour from HCC and iCCA, but we found that it also had two subtypes (CHP and CIP), indicating that the two subtypes could be used when it was difficult to distinguish H-ChC, HCC and iCCA by pathology, although larger samples are still required to verify CHP and CIP.

CHP has a metabolism-related phenotype, while CIP shows a proliferation-related phenotype

The previous results demonstrated that H-ChC is a separate biphenotypic cancer with two subtypes (CHP and CIP, Figure S4B,S4C). We identified subtype- and cluster-specific gene programs by differential expression (DE) and pathway analysis, respectively (Figure 3C-3E, Figure S4D-S4K, Tables S3,S4). Genes regulating metabolism (including the IFN-α gene) and apoptosis-related genes (e.g., the IFI6 gene) were enriched in CHP. IFN-α inhibits HIF1a signalling by suppressing FosB transcription in HCC cells, leading to a reduced glucose consumption capacity; this process establishes a high-glucose microenvironment and promotes the transcription of the T-cell costimulatory molecule CD27 in infiltrating CD8⁺ T cells through mTOR-FoxM1 signalling (30). It has been suggested that IFI6 negatively regulates endogenous apoptotic signalling pathways and TNFSF10-induced apoptosis (31-33) and exerts antiviral activity against hepatitis C virus by inhibiting EGFR signalling (34).

In contrast, CIP exhibits a gene expression pattern of epithelial cell proliferation, including expression of the tumour-associated transcription factors (TFs) SOX9 (35) and ELF3. SOX9, as a regulatory factor of epithelial progenitor cell proliferation and differentiation (36,37), participates in the branching morphogenesis of the lung epithelium by balancing proliferation and differentiation (38). ELF3 plays an important role in epithelial cell differentiation and tumorigenesis (39,40). In addition, ELF3 promotes tumour growth by activating the NF-κB pathway (41), regulates the apoptosis of intestinal epithelial cells in ulcerative colitis, and regulates the cell cycle and cell proliferation in non-small cell lung cancer and chemical lung injury (42).

To further illustrate the potential functional phenotypes of CIP and CHP, we also used the ARACNe algorithm (43) in VIPER to construct gene regulatory networks and infer TF targets based on correlation to the TF protein abundance. TFs with an adjusted P value less than 0.05 were considered significant. There were 25 significantly activated TFs in CHP: EHF, ZNF292, SOX9, JUNB, SOX4, ZNF165, ZNF281, HNF1B, EPAS1, CSRNP1, MYRF, MEF2D, FOSB, NR3C1, YY1, GTF2I, NPAS2, SMAD3, YBX3, HMG20B, ELF1, ELF3, KLF3, KLF7, KLF5. And 3 significantly activated TFs in CIP: ATF6, PA2G4, GTF3A (Figure S5D). We found that the significantly activated TFs in CHP targeted genes in metabolic pathways, whereas the significantly activated TFs in CIP targeted genes in cell proliferation (Figure S5E). These results indirectly demonstrated the potential functions of the two molecular subtypes. TFs were highly activated in the ARACNe immune-suppressive biological interaction network of CHP and CIP, which suggests that they have a certain degree of metabolic and cell proliferation effect.

A specific subpopulation enriched in CHP has characteristics of activation status

Subsequently, phenotypes across multiple tumours were analysed to determine whether there is a status that could explain the general aggressiveness of tumour cells. We identified a subpopulation named nonquiescent malignant epithelial cells in the cell cycle (NCM-Epi) as having a nonquiescent cell cycle phenotype (Figure 3E-3H and Figure S4I-S4K). Unsupervised clustering of tumour cells confirmed 9 cell clusters, and NCM-Epi was cluster 6 (including 395 cells), which was significantly enriched in the CHP subtype of H-ChC (Figure S4A-S4C and Table S4) and spanned a range of H-ChC patients (Figure S3B).

Further functional analyses were used to assess whether this population contributes to H-ChC tumorigenesis. NCM-Epi were enriched in E2F_TARGETS and G2/M_CHECKPOINT within HALLMARK terms (Figure S4K), which suggests that most of the gene signatures from cells in NCM-Epi are important factors in the cell cycle and might be crucial for cancer progression (Figure 3E). Analysis of the cell cycle revealed that the remaining clusters were in the quiescent state, while most of the cells in NCM-Epi were in the G2/M phase (with a small number being in the G1/S phase) (Figure S4I,S4J). Importantly, cells in NCM-Epi coexisted in almost all patients in this study and had elevated expression of proliferation markers, epithelial mesenchymal transition (EMT) markers and cell cycle marker genes (Figure 3F and Figure S5F).

In order to further demonstrate the immunosuppressive function of NCM-epi, TF activity for each cell was inferred using the VIPER package (43). Activity scores for NCM-epi and other malignant tumour cells were compared using Student’s t-test and the P values were adjusted using the Benjamini-Hochberg method. It was found that five TFs, DNMT1, YXB1, GTF2I, HMGA1, and PA2GA, were significantly activated in NCM-epi, which had been implicated in various biological processes and diseases (Figure 3G). Moreover, PA2GA and GTF2I have targeting relationships with genes in the PD1 signalling pathway (PTEN, JAK1, EML4, CSNK2B), with MI scores greater than 0.5 (Figure 3H). PA2GA and GTF2I were highly activated in the ARACNe immune-suppressive biological interaction network of NCM-epi, which suggests that NCM-epi has a certain degree of immunosuppressive effect and is different from proliferating cells of other types.

Collectively, these findings suggest that NCM-Epi is a tumour cell subpopulation with an immunosuppressive phenotype that is selectively increased in H-ChC. Importantly, because of their hyperactive status, NCM-Epi cells were not expected to serve as tumour stem cells. To this end, we conducted evolutionary analysis according to Slingshot (17). We found that NCM-Epi could act as the intermediate status cell subgroup (Figure S4F,S4G), which is inconsistent with the dedifferentiated phenotype of cancer stem cells with promiscuous signals (44,45). However, subsequent work with an expanded cohort size is needed to validate these results.

H-ChC shows greater exhausted CD8⁺ T-cell dysfunction than HCC and iCCA

To describe the immune cell status in the H-ChC TIME, immune cells from 7 H-ChC samples were extracted in the cohort (N=28,760 cells), with immune cells from HCC samples (N=20,235 cells) and iCCA samples (N=18,437 cells) used as controls (Figure 4A and Figure S6A). Taking into account the heterogeneity in cellular composition among patients, we estimated cell proportions using a bias-corrected and accelerated bootstrap algorithm (BCa) (16), which uses bootstrap resampling (1,000 times), to calculate the standard error associated with estimates of cell type proportion and conduct statistical analyses of the three tumour types (Figure 4B).

Figure 4 Greater exhausted CD8⁺ T-cell dysfunction was detected in H-ChC. (A) t-SNE projections of the nonepithelial subsets: B-cell compartment, T-cell/NK compartment myeloid compartment and stromal compartment of H-ChC, HCC and iCCA. All cells are coloured according to their cell type. (B) Heatmap showing the proportions of nonepithelial cell subsets in different tumours. Colour indicates the scaled cellular proportion. (C) The proportions of selected nonepithelial cell subsets. The percentage is presented along the y-axis (bootstrap), and the different tumour types are presented along the x-axis. The shaded areas represent the upper and lower quantile bootstrap cell proportions. (D) The incoming and outgoing interaction strengths of different cell types in H-ChC (right) and iCCA (left). In contrast to the findings for iCCA, CD8⁺ Tex cells become one of the main sources and targets of H-ChC. (E) The ranked differential CD8⁺ Tex cells with other cell crosstalk (H-ChC versus iCCA). H-ChC, hepatocholangiocarcinoma; HCC, hepatocellular carcinoma; iCCA, intrahepatic cholangiocarcinoma; Tex, exhausted T cell; EC, endothelial cell; t-SNE, T-distribution stochastic neighbour embedding.

The immune cell abundance can be divided into three patterns: H-ChC enriched, HCC enriched, and iCCA enriched. CD8⁺ Tex, CD4⁺ follicular helper T cell (Tfh), CD4⁺ Treg, CD8⁺ MAIT, CD8⁺ prolifer ative T cell (Tprf), and CD8⁺ tissue-resident memory T cell (Trm) cells were enriched in H-ChC samples, with CXCL10⁺ macrophages, inflammatory carcinoma-associated fibroblasts (iCAFs), lymphatic ECs, and SPP1⁺ macrophages being mostly enriched in iCCA samples (Figure 4B,4C, Figure S6B). Native T cell (Tn) cells, NK cells, Tip ECs, and SLC40A1⁺ C1QC⁺ macrophages were enriched in HCC, consistent with previously reported results (46,47). These results suggest a tumour-specific distribution of certain immune cells. Interestingly, these cells were present simultaneously in the TIME, suggesting a common tumour education role of immune cells (Figure S6A and Figure 4B). Immunosuppressive cells accounted for a substantial part of H-ChC, comprising the surge of CD8⁺ Tex cells and neutrophils with potential tumour-promoting effects, which was similar to a previous study (Figure 4C and Figure S6B) (48).

Given the three patterns of immune cell abundance, we hypothesized that different tumour-specific distributions may function distinctly in the three tumours. We performed a ligand-receptor crosstalk analysis and compared the three tumour types (Figure 4D and Figure S6C). The total number and strength of interactions were highest in H-ChC (Figure S6D), followed by iCCA (weakest interaction strength) and HCC (lowest number of interactions). Notably, the signal fluctuation (interaction strength of incoming and outgoing) of CD8⁺ Tex cells was the strongest, which was obviously higher than that of HCC and iCCA cells (Figure 4D and Figure S6C). CXCL10⁺ macrophages and lymphatic ECs ranked high among the crosstalk cell types for CD8⁺ Tex cells in the H-ChC samples (Figure 4E and Figure S7A-S7D). CXCL10⁺ macrophages highly express CXCR4 (49) and CCR1, key molecules that promote tumour stem cell-like properties and influence stem cell proliferation (50,51) (Figure S7A,S7B, Table S5). CCR5, a chemotactic receptor, participates in the migration of T lymphocytes to the site of infection (52). In contrast to CXCL10⁺ macrophages, we found that lymphatic ECs exhibited higher ACKR2 expression than CXCR4 or CCR1 expression. By removing chemokines from tissues, lymphatic vessel surfaces, and the placenta, ACKR2 plays an important role in the resolution (termination) of the inflammatory response and in the regulation of the adaptive immune response (53). These data suggested that the immunological inhibition phenotype of CXCL10⁺ macrophages may be shaped by the suppressive antitumour immune response in the TIME (Figure 4E).

scRNA-seq detected GZMK⁺ CD8⁺ Tex cells in H-ChC, HCC and iCCA

Previous studies have demonstrated that PD-1 pathway blockade can terminally reduce the dysfunctional state of CD8⁺ Tex cells, suggesting an extensive mechanism for immune checkpoint blockade (ICB) therapy (54-56). However, it is still unknown how most T cells are deprived of their tumoricidal action and become CD8⁺ Tex cells, which play a major regulatory role in H-ChC. To gain further insight into the molecular heterogeneity of CD8⁺ Tex cells in response to ICB therapy, we clustered 6,852 CD8⁺ Tex cells into two broad cell states, GZMK⁺ CD8⁺ Tex cells and terminal CD8⁺ Tex cells (Figure S5A-S5C). GZMK⁺ CD8⁺ Tex cells were characterized by high expression of GZMK, while terminal CD8⁺ Tex cells were identified by low expression of GZMK, which was an effector signature (Figure S8A-S8C). To address whether tumour regression was mainly driven by these two CD8⁺ Tex cell states (Figure 5A), we observed a higher infiltration of terminal CD8⁺ Tex cells but a higher strength of interactions for GZMK⁺ CD8⁺ Tex cells in the TIME of H-ChC (Figure 5B-5D).

Figure 5 Cell communication network of GZMK⁺ CD8⁺ Tex cells and terminal Tex cells. (A) Expression levels of selected genes in CD8⁺ Tex subpopulations in H-ChC tumours. *, P<0.05; **, P<0.01; ***, P<0.001. (B) Interaction strengths between selected cell types and CD8⁺ Tex subpopulations in H-ChC tumours. (C) Interaction strengths between selected cell types and CD8⁺ Tex subpopulations in HCC tumours. (D) Interaction strengths between selected cell types and CD8⁺ Tex subpopulations in iCCA tumours. H-ChC, hepatocholangiocarcinoma; HCC, hepatocellular carcinoma; iCCA, intrahepatic cholangiocarcinoma; Tex, exhausted T cell.

Since multiple mouse studies have demonstrated that anti-PD-1 could also target other types of Tex cells to enhance antitumour immunity (54-56). We hypothesized that GZMK⁺ CD8⁺ Tex cells might interact with tumour cells to preferentially decrease progression. We detected the different regulatory networks of GZMK⁺ CD8⁺ Tex cells based on CHP, CIP and NCM-Epi (Figure S8D,S8E). Notably, the interaction between GZMK⁺ CD8⁺ Tex cells and malignant cells in the TIME of the three tumour types was strongest for NCM-Epi (Figure S8E).

A deeper inspection of specific cell-ligand pairs that were enriched in CD8⁺ GZMK⁺ Tex cells and NCM-Epi was conducted. Using partial correlation coefficients (PCCs) and the gene set enrichment analysis (GSEA) algorithm (57), we identified immune-related gene sets in NCM-Epi. All pairs were ranked based on RS scores and then subjected to enrichment analysis, in which COL4A2-CD44 is a specific cell ligand pair that has been found to be enriched in CD8⁺ GZMK⁺ Tex cells and NCM-Epi. COL4A2 is a collagen protein that plays a role in extracellular matrix assembly and remodelling, and has been implicated in promoting tumour invasiveness and angiogenesis. CD44 is a transmembrane glycoprotein that is involved in cell-cell interactions and signalling, and has been shown to be upregulated in a variety of cancer types, where it promotes tumour cell survival, migration, invasion and metastasis. The interaction between COL4A2 and CD44 may thus contribute to the infiltration and control of immune cells in the tumour microenvironment, and potentially influence the prognosis of cancer patients. Together, our findings suggest that these GZMK⁺ CD8⁺ Tex cells may exist in an immunosuppressed TIME characterized by the NCM-Epi subpopulation.

GZMK⁺ CD8⁺ Tex cells may play a key role in clinical immunotherapy

Given that most tumours regress (albeit incompletely) during PD-1 immunotherapy, it is extremely challenging to obtain serial clinical tumour biopsy samples from a single patient with advanced H-ChC cancer before and during treatment (58,59). However, we obtained a sample before the treatment time point from the tumour of a patient (S16) who received targeted treatment and anti-PD-1 immunotherapy (Figure S9A). Therefore, we used scRNA-seq to investigate the temporal dynamics of immune cells and to characterize the features associated with immunotherapy.

We obtained the single-cell transcriptome data of 33,918 immune cells, including myeloid, mast, T and B cells, through high-quality screening (Figure S9B,S9C, Table S6). We also identified cell subtypes for all major cell types (Figure 6A). The abundances of major types of immune cells changed dynamically after treatment: the numbers of lymphocytes increased, those of myeloid cells decreased (Figure S9C), and the tumour cells no longer existed after treatment.

Figure 6 CD8⁺ Tex cells may act as a potential driver for H-ChC development. (A) t-SNE plot of all cells coloured according to cell type. (B) The cellular compositions of T-cell subtypes in H-ChC tumours pre- and posttreatment (left). Box plots showing the alterations in the effector memory score of CD8⁺ Tex cells (right, STAR Methods). (C) Expression levels of selected genes in CD8⁺ Tex subpopulations in H-ChC tumours pre- and posttreatment. (D) The difference between the incoming signal strength of H-ChC pre- and posttreatment. (E) The cellular compositions of myeloid (left) and immune (right) cell subtypes in H-ChC tumours pre- and posttreatment. (F) Kaplan-Meier analysis of follow-up time and PFS in an independent cohort of H-ChC patients. ***, P<0.001. TIME, tumour immune microenvironment; H-ChC, hepatocholangiocarcinoma; HCC, hepatocellular carcinoma; iCCA, intrahepatic cholangiocarcinoma; Tex, exhausted T cell; t-SNE, T-distribution stochastic neighbour embedding; PFS, progression-free survival.

Considering the significance of T cells in anti-programmed cell death ligand 1 (anti-PD-L1) therapy, higher levels of T cells were observed after immunotherapy (Figure 6B and Figure S9C). Subsequently, we investigated the temporal transcriptional dynamics of the different T-cell subpopulations induced by treatment (Figure 6C). Transcriptome analysis showed that the expression of effector genes such as IFNG, GZMH, and GZMB was upregulated after treatment, whereas the expression of exhaustion-related genes in CD8⁺ Tex cells, including CTLA4 and TIGIT, was downregulated (Figure 6B,6C; Table S7). The phenotypic transformation in CD8⁺ Tex cells was also manifested by a decreased exhaustion degree and an increased effector memory score after treatment (Figure 6B-6D). In addition, H-ChC showed high infiltration of GZMK⁺ CD8⁺ Tex cells after treatment (Figure 6E), considering that PD-1 blockade may preferentially increase their proportions and further expand tumour-reactive T cells following treatment (Figure 6E and Figure S9D), which may indicate that immunotherapy would benefit patients with advanced H-ChC.

In our analysis, we found that the deficiency of tumour cells after immunotherapy was just one feature of H-ChC, with stronger levels of interaction between GZMK⁺ CD8⁺ Tex cells and other proinflammatory cells after immunotherapy (Figure S9E) (60,61). Thus, we recruited 10 patients with advanced H-ChC into one of our real-world studies, where they were treated with PD-1 and lenvatinib. With a median follow-up time of 23.3 months, the median progression-free survival (PFS) was 10.1 months (IQR, 3.0–17.2 months) (Figure 6F). Regarding the best overall response, according to RECIST version 1.1, 6 patients achieved PR, 3 patients achieved SD, and 1 achieved PD. The ORR was 60%, and the DCR was 90%, which suggests that H-ChC is more sensitive than HCC or iCCA to immunotherapy. In further studies, H-ChC should be considered a separate tumour and observed for its response after immunotherapy.

Discussion

H-ChC is a rare primary hepatic carcinoma with varying degrees of differentiation of hepatic carcinoma and cholangiocarcinoma due to its substantial heterogeneity in morphological, immunophenotypic, molecular, radiological, and clinical features; this heterogeneity is a challenge to overcome. Our previous work, a CNV analysis based on exome data, revealed substantial intratumour heterogeneity within H-ChC and suggested that H-ChC is monoclonal (4,5). The factors that determine HCC and iCCA phenotypes in H-ChC remain to be explored (62-64).

Elucidation of the biological complexity that cannot be described by bulk-level subtyping was sought here, and the findings support the view that H-ChC is an independent biphenotypic cancer, and the two components we identified, CHP and CIP, had similar functions in iCCA and HCC, respectively. However, CHP highly expressed iCCA epithelial markers, and CIP highly expressed HCC markers (Figure 3B). ScRNA-seq affords a unique opportunity to assess the potential for plasticity between subtypes, and these subtypes may not have been captured by previous techniques (22,23,65-69), particularly the specific malignant cell state, NCM-Epi. The tumour stem cell populations may correspond to the upstream of NCM-Epi and result in the loss of the undifferentiated state of stem cells and/or hepatocytes or cholangiocyte differentiation and then generate multiple epithelial morphologies (e.g., hepatocytes and/or cholangiocytes) owing to the promiscuous signalling property.

In our study, we also characterized the immunogenicity of biphenotypic carcinoma, which may drive future work to ultimately decode the different formations of TIME properties in biphenotypic carcinoma. Nevertheless, the consistent presence of cold immune phenotypes (70) with latent impacts on the determinants of the dual phenotype of H-ChC was distinguished. We obtained pre- and posttreatment tumour samples from a patient (P017) receiving anti-PD-1 immunotherapy and found that the CD8⁺ Tex cells from this patient exhibited phenotypic alterations after therapy, with a decreased exhaustion score and increased effector memory score. These results may indicate that many antigens recognized by the immune system may be present on tumour cells in biphenotypic carcinoma, further suggesting enhanced effector properties after PD-L1 blockade treatment.

Furthermore, CD8⁺ Tex cells were significantly enriched in H-ChC, followed by iCCA and HCC. The patient had increased levels of GZMK⁺ CD8⁺ Tex cells after PD-1 treatment, and we posit that inhibiting the interaction between NCM-Epi and GZMK⁺ CD8⁺ Tex cells is a critical step towards improving the therapeutic treatment of H-ChC. In addition, the higher expression of TIGIT and CTAL4 in H-ChC than in HCC and iCCA implies that immunotherapy is more suitable for biphenotypic carcinoma. Moreover, potential treatment methods for biphenotypic carcinoma should be discussed after collecting tumours from new chemotherapy or immunotherapy groups, rather than simply applying treatments for HCC and iCCA, although both HCC and iCCA treatments may be effective.

Mechanistic studies employing lineage tracing models can help locate the origins of cells of H-ChC populations. However, our data describe the landscape of H-ChC, HCC, and iCCA tumours in detail and reveal the impacts of the ecological niche and malignant tumours on distinct subtypes, thus offering valuable resources for future studies of specific tumours. Such a detailed understanding will pave the way for the development of targeted immunotherapeutic strategies for specific components of TIMEs.

Conclusions

Here, we used scRNA-seq data analysis to systematically compare the transcriptomic signatures of different cell populations in HCC, iCCA and H-ChC, which revealed the molecular mechanisms underlying the coexistence of HCC and iCCA components in H-ChC.

Acknowledgments

We thank the Department of Liver Surgery at Peking Union Medical College Hospital for providing clinical tumour biopsy samples collected from patients with H-ChC before and during treatment, which are extremely challenging to obtain. And Beijing SeekGene BioSciences Co., Ltd. (Beijing, China) for assistance in scRNA-seq related experiments and tissue dissociation.

Funding: This work was supported by the CAMS Innovation Fund for Medical Sciences (CIFMS) (2021-I2M-1-061 and 2021-1-I2M-003), the CSCO-hengrui Cancer Research Fund (Y-HR2019-0239), the CSCO-MSD Cancer Research Fund (Y-MSDZD2021-0213), the National Ten-thousand Talent Program (to J.L.), the Young Scientists Fund of the National Natural Science Foundation of China (82303720), the Beijing Natural Science Foundation (7234381), and the Fundamental Research Funds for the Central Universities (3332023011 to J.L.).

Footnote

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://hbsn.amegroups.com/article/view/10.21037/hbsn-23-400/coif). H.Z. serves as an unpaid Deputy Editor-in-Chief of Hepatobiliary Surgery and Nutrition, and he reports that this work was supported by the CAMS Innovation Fund for Medical Sciences (CIFMS) (2021-I2M-1-061 and 2021-1-I2M-003), the CSCO-hengrui Cancer Research Fund (Y HR2019-0239), the CSCO-MSD Cancer Research Fund (Y-MSDZD2021-0213) and the National Ten-thousand Talent Program. The other authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The studies involving human participants were reviewed and approved by the institutional review board and ethics committee at Peking Union Medical College Hospital (No. JS-1391). The patients/participants provided their written informed consent to participate in this study.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Bray F, Ferlay J, Soerjomataram I, et al. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2018;68:394-424. [Crossref] [PubMed]
2. Sciarra A, Park YN, Sempoux C. Updates in the diagnosis of combined hepatocellular-cholangiocarcinoma. Hum Pathol 2020;96:48-55. [Crossref] [PubMed]
3. Lee JH, Chung GE, Yu SJ, et al. Long-term prognosis of combined hepatocellular and cholangiocarcinoma after curative resection comparison with hepatocellular carcinoma and cholangiocarcinoma. J Clin Gastroenterol 2011;45:69-75. [Crossref] [PubMed]
4. Wang A, Wu L, Lin J, et al. Whole-exome sequencing reveals the origin and evolution of hepato-cholangiocarcinoma. Nat Commun 2018;9:894. [Crossref] [PubMed]
5. Xue R, Chen L, Zhang C, et al. Genomic and Transcriptomic Profiling of Combined Hepatocellular and Intrahepatic Cholangiocarcinoma Reveals Distinct Molecular Subtypes. Cancer Cell 2019;35:932-947.e8. [Crossref] [PubMed]
6. Butler A, Hoffman P, Smibert P, et al. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol 2018;36:411-20. [Crossref] [PubMed]
7. Aliverti E, Tilson JL, Filer DL, et al. Projected t-SNE for batch correction. Bioinformatics 2020;36:3522-7. [Crossref] [PubMed]
8. Li Q, Wang R, Yang Z, et al. Molecular profiling of human non-small cell lung cancer by single-cell RNA-seq. Genome Med 2022;14:87. [Crossref] [PubMed]
9. Gao R, Bai S, Henderson YC, et al. Delineating copy number and clonal substructure in human tumors from single-cell transcriptomes. Nat Biotechnol 2021;39:599-608. [Crossref] [PubMed]
10. Sharma G, Colantuoni C, Goff LA, et al. projectR: an R/Bioconductor package for transfer learning via PCA, NMF, correlation and clustering. Bioinformatics 2020;36:3592-3. [Crossref] [PubMed]
11. Gonzalez H, Mei W, Robles I, et al. Cellular architecture of human brain metastases. Cell 2022;185:729-745.e20. [Crossref] [PubMed]
12. Chan JM, Quintanal-Villalonga Á, Gao VR, et al. Signatures of plasticity, metastasis, and immunosuppression in an atlas of human small cell lung cancer. Cancer Cell 2021;39:1479-1496.e18. [Crossref] [PubMed]
13. Yang C, Zhang S, Cheng Z, et al. Multi-region sequencing with spatial information enables accurate heterogeneity estimation and risk stratification in liver cancer. Genome Med 2022;14:142. [Crossref] [PubMed]
14. Yu G, Wang LG, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012;16:284-7. [Crossref] [PubMed]
15. Ma L, Hernandez MO, Zhao Y, et al. Tumor Cell Biodiversity Drives Microenvironmental Reprogramming in Liver Cancer. Cancer Cell 2019;36:418-430.e6. [Crossref] [PubMed]
16. Cao Y, Lin Y, Ormerod JT, et al. scDC: single cell differential composition analysis. BMC Bioinformatics 2019;20:721. [Crossref] [PubMed]
17. Street K, Risso D, Fletcher RB, et al. Slingshot: cell lineage and pseudotime inference for single-cell transcriptomics. BMC Genomics 2018;19:477. [Crossref] [PubMed]
18. WangYWangRZhangSiTALK: an R Package to Characterize and Illustrate Intercellular Communication.bioRxiv 2019. doi: https://doi.org/.
19. Jin S, Guerrero-Juarez CF, Zhang L, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun 2021;12:1088. [Crossref] [PubMed]
20. Cheng S, Li Z, Gao R, et al. A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells. Cell 2021;184:792-809.e23. [Crossref] [PubMed]
21. Zhang Q, He Y, Luo N, et al. Landscape and Dynamics of Single Immune Cells in Hepatocellular Carcinoma. Cell 2019;179:829-845.e20. [Crossref] [PubMed]
22. Zhang X, Peng L, Luo Y, et al. Dissecting esophageal squamous-cell carcinoma ecosystem by single-cell transcriptomic analysis. Nat Commun 2021;12:5291. [Crossref] [PubMed]
23. Zhang Y, Chen H, Mo H, et al. Single-cell analyses reveal key immune cell subsets associated with response to PD-L1 blockade in triple-negative breast cancer. Cancer Cell 2021;39:1578-1593.e8. [Crossref] [PubMed]
24. Liu B, Li C, Li Z, et al. An entropy-based metric for assessing the purity of single cell populations. Nat Commun 2020;11:3155. [Crossref] [PubMed]
25. Liu Z, Yu S, Ye S, et al. Keratin 17 activates AKT signalling and induces epithelial-mesenchymal transition in oesophageal squamous cell carcinoma. J Proteomics 2020;211:103557. [Crossref] [PubMed]
26. Teschendorff AE, Enver T. Single-cell entropy for accurate estimation of differentiation potency from a cell's transcriptome. Nat Commun 2017;8:15599. [Crossref] [PubMed]
27. Xue R, Zhang Q, Cao Q, et al. Liver tumour immune microenvironment subtypes and neutrophil heterogeneity. Nature 2022;612:141-7. [Crossref] [PubMed]
28. Ma L, Wang L, Khatib SA, et al. Single-cell atlas of tumor cell evolution in response to therapy in hepatocellular carcinoma and intrahepatic cholangiocarcinoma. J Hepatol 2021;75:1397-408. [Crossref] [PubMed]
29. Pelka K, Hofree M, Chen JH, et al. Spatially organized multicellular immune hubs in human colorectal cancer. Cell 2021;184:4734-4752.e20. [Crossref] [PubMed]
30. Hu B, Yu M, Ma X, et al. IFNα Potentiates Anti-PD-1 Efficacy by Remodeling Glucose Metabolism in the Hepatocellular Carcinoma Microenvironment. Cancer Discov 2022;12:1718-41. [Crossref] [PubMed]
31. Cheriyath V, Glaser KB, Waring JF, et al. G1P3, an IFN-induced survival factor, antagonizes TRAIL-induced apoptosis in human myeloma cells. J Clin Invest 2007;117:3107-17. [Crossref] [PubMed]
32. Qi Y, Li Y, Zhang Y, et al. IFI6 Inhibits Apoptosis via Mitochondrial-Dependent Pathway in Dengue Virus 2 Infected Vascular Endothelial Cells. PLoS One 2015;10:e0132743. [Crossref] [PubMed]
33. Tahara E Jr, Tahara H, Kanno M, et al. G1P3, an interferon inducible gene 6-16, is expressed in gastric cancers and inhibits mitochondrial-mediated apoptosis in gastric cancer cell line TMK-1 cell. Cancer Immunol Immunother 2005;54:729-40. [Crossref] [PubMed]
34. Meyer K, Kwon YC, Liu S, et al. Interferon-α inducible protein 6 impairs EGFR activation by CD81 and inhibits hepatitis C virus infection. Sci Rep 2015;5:9012. [Crossref] [PubMed]
35. Angelici B, Shen L, Schreiber J, et al. An AAV gene therapy computes over multiple cellular inputs to enable precise targeting of multifocal hepatocellular carcinoma in mice. Sci Transl Med 2021;13:eabh4456. [Crossref] [PubMed]
36. van Gastel N, Stegen S, Eelen G, et al. Lipid availability determines fate of skeletal progenitor cells via SOX9. Nature 2020;579:111-7. [Crossref] [PubMed]
37. Villani V, Thornton ME, Zook HN, et al. SOX9+/PTF1A+ Cells Define the Tip Progenitor Cells of the Human Fetal Pancreas of the Second Trimester. Stem Cells Transl Med 2019;8:1249-64. [Crossref] [PubMed]
38. Li L, Feng J, Zhao S, et al. SOX9 inactivation affects the proliferation and differentiation of human lung organoids. Stem Cell Res Ther 2021;12:343. [Crossref] [PubMed]
39. Chen H, Chen W, Zhang X, et al. E26 transformation (ETS)-specific related transcription factor-3 (ELF3) orchestrates a positive feedback loop that constitutively activates the MAPK/Erk pathway to drive thyroid cancer. Oncol Rep 2019;41:570-8. [Crossref] [PubMed]
40. Yachida S, Wood LD, Suzuki M, et al. Genomic Sequencing Identifies ELF3 as a Driver of Ampullary Carcinoma. Cancer Cell 2016;29:229-40. [Crossref] [PubMed]
41. Li L, Miao X, Ni R, et al. Epithelial-specific ETS-1 (ESE1/ELF3) regulates apoptosis of intestinal epithelial cells in ulcerative colitis via accelerating NF-κB activation. Immunol Res 2015;62:198-212. [Crossref] [PubMed]
42. He D, Wang D, Lu P, et al. Single-cell RNA sequencing reveals heterogeneous tumor and immune cell populations in early-stage lung adenocarcinomas harboring EGFR mutations. Oncogene 2021;40:355-68. [Crossref] [PubMed]
43. Huang C, Chen L, Savage SR, et al. Proteogenomic insights into the biology and treatment of HPV-negative head and neck squamous cell carcinoma. Cancer Cell 2021;39:361-379.e16. [Crossref] [PubMed]
44. Zhang C, Li J, Cheng Y, et al. Single-cell RNA sequencing reveals intrahepatic and peripheral immune characteristics related to disease phases in HBV-infected patients. Gut 2023;72:153-67. [Crossref] [PubMed]
45. Zhang Z, Wang ZX, Chen YX, et al. Integrated analysis of single-cell and bulk RNA sequencing data reveals a pan-cancer stemness signature predicting immunotherapy response. Genome Med 2022;14:45. [Crossref] [PubMed]
46. Duan M, Goswami S, Shi JY, et al. Activated and Exhausted MAIT Cells Foster Disease Progression and Indicate Poor Outcome in Hepatocellular Carcinoma. Clin Cancer Res 2019;25:3304-16. [Crossref] [PubMed]
47. Song G, Shi Y, Zhang M, et al. Global immune characterization of HBV/HCV-related hepatocellular carcinoma identifies macrophage and T-cell subsets associated with disease progression. Cell Discov 2020;6:90. [Crossref] [PubMed]
48. Zilionis R, Engblom C, Pfirschke C, et al. Single-Cell Transcriptomics of Human and Mouse Lung Cancers Reveals Conserved Myeloid Populations across Individuals and Species. Immunity 2019;50:1317-1334.e10. [Crossref] [PubMed]
49. Fan QM, Jing YY, Yu GF, et al. Tumor-associated macrophages promote cancer stem cell-like properties via transforming growth factor-beta1-induced epithelial-mesenchymal transition in hepatocellular carcinoma. Cancer Lett 2014;352:160-8. [Crossref] [PubMed]
50. Meng J, Wang L, Hou J, et al. CCL23 suppresses liver cancer progression through the CCR1/AKT/ESR1 feedback loop. Cancer Sci 2021;112:3099-110. [Crossref] [PubMed]
51. Su S, Mukaida N, Wang J, et al. Inhibition of immature erythroid progenitor cell proliferation by macrophage inflammatory protein-1alpha by interacting mainly with a C-C chemokine receptor, CCR1. Blood 1997;90:605-11.
52. Sharapova TN, Romanova EA, Sashchenko LP, et al. Tag7-Mts1 Complex Induces Lymphocytes Migration via CCR5 and CXCR3 Receptors. Acta Naturae 2018;10:115-20.
53. Mantovani A. Wandering pathways in the regulation of innate immunity and inflammation. J Autoimmun 2017;85:1-5. [Crossref] [PubMed]
54. Im SJ, Hashimoto M, Gerner MY, et al. Defining CD8+ T cells that provide the proliferative burst after PD-1 therapy. Nature 2016;537:417-21. [Crossref] [PubMed]
55. Miller BC, Sen DR, Al Abosy R, et al. Subsets of exhausted CD8(+) T cells differentially mediate tumor control and respond to checkpoint blockade. Nat Immunol 2019;20:326-36. [Crossref] [PubMed]
56. Siddiqui I, Schaeuble K, Chennupati V, et al. Intratumoral Tcf1(+)PD-1(+)CD8(+) T Cells with Stem-like Properties Promote Tumor Control in Response to Vaccination and Checkpoint Blockade Immunotherapy. Immunity 2019;50:195-211.e10. [Crossref] [PubMed]
57. Li Y, Jiang T, Zhou W, et al. Pan-cancer characterization of immune-related lncRNAs identifies potential oncogenic biomarkers. Nat Commun 2020;11:1000. [Crossref] [PubMed]
58. Camidge DR, Dziadziuszko R, Peters S, et al. Updated Efficacy and Safety Data and Impact of the EML4-ALK Fusion Variant on the Efficacy of Alectinib in Untreated ALK-Positive Advanced Non-Small Cell Lung Cancer in the Global Phase III ALEX Study. J Thorac Oncol 2019;14:1233-43. [Crossref] [PubMed]
59. Soria JC, Ohe Y, Vansteenkiste J, et al. Osimertinib in Untreated EGFR-Mutated Advanced Non-Small-Cell Lung Cancer. N Engl J Med 2018;378:113-25. [Crossref] [PubMed]
60. Lavin Y, Kobayashi S, Leader A, et al. Innate Immune Landscape in Early Lung Adenocarcinoma by Paired Single-Cell Analyses. Cell 2017;169:750-765.e17. [Crossref] [PubMed]
61. Zhang L, Li Z, Skrzypczynska KM, et al. Single-Cell Analyses Inform Mechanisms of Myeloid-Targeted Therapies in Colon Cancer. Cell 2020;181:442-459.e29. [Crossref] [PubMed]
62. Gerlinger M, Rowan AJ, Horswell S, et al. Intratumor heterogeneity and branched evolution revealed by multiregion sequencing. N Engl J Med 2012;366:883-92. [Crossref] [PubMed]
63. Miao R, Luo H, Zhou H, et al. Identification of prognostic biomarkers in hepatitis B virus-related hepatocellular carcinoma and stratification by integrative multi-omics analysis. J Hepatol 2014;61:840-9. [Crossref] [PubMed]
64. Xue R, Li R, Guo H, et al. Variable Intra-Tumor Genomic Heterogeneity of Multiple Lesions in Patients With Hepatocellular Carcinoma. Gastroenterology 2016;150:998-1008. [Crossref] [PubMed]
65. Guo L, Yi X, Chen L, et al. Single-Cell DNA Sequencing Reveals Punctuated and Gradual Clonal Evolution in Hepatocellular Carcinoma. Gastroenterology 2022;162:238-52. [Crossref] [PubMed]
66. Liu B, Hu X, Feng K, et al. Temporal single-cell tracing reveals clonal revival and expansion of precursor exhausted T cells during anti-PD-1 therapy in lung cancer. Nat Cancer 2022;3:108-21. [Crossref] [PubMed]
67. Liu Y, Zhang Q, Xing B, et al. Immune phenotypic linkage between colorectal cancer and liver metastasis. Cancer Cell 2022;40:424-437.e5. [Crossref] [PubMed]
68. Wang D, Hou S, Zhang L, et al. iMAP: integration of multiple single-cell datasets by adversarial paired transfer networks. Genome Biol 2021;22:63. [Crossref] [PubMed]
69. Wang X, He Y, Zhang Q, et al. Direct Comparative Analyses of 10X Genomics Chromium and Smart-seq2. Genomics Proteomics Bioinformatics 2021;19:253-66. [Crossref] [PubMed]
70. Liu Y, Yan X, Zhang F, et al. TCR-T Immunotherapy: The Challenges and Solutions. Front Oncol 2021;11:794183. [Crossref] [PubMed]