Tracking the important role of JUNB in hepatocellular carcinoma by single-cell sequencing analysis.

Hepatocellular carcinoma (HCC) is the most commonly diagnosed liver cancer, accounting for ~90% of all primary malignancy of the liver. Although various medical treatments have been used as systemic therapies, patient survival time may be extended by only a few months. Moreover, the underlying mechanisms of HCC development and progression remain poorly understood. In the present study, the single-cell transcriptome of one in vivo HCC tumor sample, two in vitro HCC cell lines and normal peripheral blood mononuclear cells were analysed in order to identify the potential mechanism underlying the development and progression of HCC. Interestingly, JunB proto-oncogene was identified to serve a role in the immune response and in development and progression of HCC, potentially contributing to the development of novel therapeutics for HCC patients. Copyright: © Yan et al.

Introduction

Hepatocellular carcinoma (HCC) is a malignant tumor with high incidence and the number of novel cases is expected to increase by 1 million every year in the next decade (1). The main causative factors include alcoholic liver disease, hepatitis B and hepatitis C viral infections, and non-alcoholic fatty liver disease (2). However, the pathogenesis and exact molecular mechanism of HCC are not fully understood. Although international guidelines suggest HCC screening for patients with cirrhosis, regular monitoring presents various limitations in clinical practice (3). The large number of undiagnosed patients leads to a low HCC monitoring rate and high late diagnosis rate, causing the patients with HCC to be diagnosed only when the tumor exhibits a large size, leading to poor prognosis (4).

During tumor development and progression, multiple cell types interact with tumor cells, including astrocytes, B cells, lymphocytes, macrophages, monocytes, natural killer cells and T cells, constituting the tumor immune microenvironment (5,6). These cells, together with the fibroblasts, vascular endothelial cells and other factors, which are collectively called the tumor stroma, as well as the extracellular matrix, oxygen levels and pH values, constitute the tumor microenvironment (7,8). Notably, the interactions between immune cells and tumor cells affect the growth and remodeling of the tumor microenvironment (9). Immune cells can stimulate tumor cells to secrete cytokines, which mediate the tumor growth by promoting the growth of new blood and lymphatic vessels (6). Different cell types may serve pro- and anti-tumorigenic roles in the tumor microenvironment. For example, targeting T cell activation is considered as an important novel strategy to repress tumor growth (10,11). S100A4 has been shown to be an oncogene able to promote inflammation (12) and affect angiogenesis (13). Accumulating evidence showed that expression of S100A4 in tumor cells is related to the tumor-associated T cell deficiency (14). The rich blood supply and unique sinusoid structure of the liver provide a plastic environment for the formation and function of the tumor immune microenvironment (15). Therefore, it is of great significance to study the molecular characteristics and intercellular interactions in the HCC immune microenvironment.

Massively parallel sequencing data have provided novel insights in the field of cancer research. In particular, RNA-sequencing (RNA-seq) has been used to detect genomic mutations and rearrangement signatures in the human genome and transcriptome. However, conventional bulk RNA-seq can only provide the average expression signal of transcripts in the whole tumor tissue, without considering the tumor heterogeneity. By contrast, single-cell RNA-seq may facilitate the identification of complex and rare cells populations, thus allowing investigation of the tumor immune microenvironment (16), especially in HCC. For example, using single-cell RNA-seq, a previous study identified 11 HCC-related T cell subpopulations, which provided valuable insights for the understanding of the cancer immune microenvironment (17). In addition, a previous study has described the molecular characteristics of immune cells that infiltrate HCC to determine whether certain types of drugs may be effective against liver cancer (18).

Moreover, chromatin immunoprecipitation (ChIP)-seq results and protein-protein interaction (PPI) network may facilitate the detection of gene regulatory networks and interaction events, such as the bindings between transcription factors (TFs) and promoters. The present study compared the single-cell RNA-seq data of normal peripheral blood mononuclear cells (PBMCs) with that of in vivo tumor cells and in vitro cell lines using single cell classification and identification. By integrating differential expression analysis, ChIP-seq data and PPI networks, the present results suggested that the JunB proto-oncogene (JUNB) may serve an important role in the development and progression of HCC and the immune response. In addition, apolipoprotein A2 (APOA2), which encodes a genetically susceptible protein in HCC (19), was found to exhibit the same expression pattern as JUNB. The present results may contribute to the identification of novel therapeutic targets for the treatment of HCC patients.

Materials and methods

Data collection

In vivo tumor cells were isolated from a patient who had undergone resection at the National Institutes of Health (NIH) Clinical Center. The tissue acquisition procedures were approved by the Institutional Review Board of The NIH (20). In total, two in vitro cell lines (HuH1 and HuH7) from The Health Science Research Resources Bank (cat. nos. JCRB0199 and JCRB0403) were pooled and used for 10× Genomics single-cell RNA-seq. These data were collected and downloaded from the GEO database (database no. GSE103867) (20). Single-cell data of PBMCs were downloaded from the GEO database (database no. GSE111360) (21). Gene expression profile and clinical information of a cohort of 360 patients with HCC were collected from The Cancer Genome Atlas (TCGA; http://cancergenome.nih.gov/). The PPI network was obtained using STRING (v11.0) with only highly strong interactions (score, 0.4) being used (22).

Preprocessing for 10× Genomics single-cell RNA-seq data

Seurat v2.1 (http://satijalab.org/seurat/) was used to analyze the 10× Genomics data (23). Genes whose expression was detected in ≥3 cells and cells with ≥10 genes were used in this study. Variable genes were identified using cutoffs (x.low.cutoff = 0.05; y.cutoff = 0.1). The top 20 principal components were used in the clustering analysis (resolution = 0.6). Gene expression levels were quantified using the unique molecular identifier counts. Dimensionality reduction was based on the t-SNE algorithm. Subsequently, cell populations were clustered by principal component analysis.

Differential expression gene and pathway enrichment analysis

Log2Fold-Change represented the ratio of gene expression between one cluster of cells and all the other cells. P-values were calculated using the negative binomial test and adjusted by the Benjamini-Hochberg method. Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using DAVID (version 6.8) (24).

Correlation analysis

Gene expression of 360 patients with HCC, obtained from the TCGA database, was used as an independent external test set to validate the putative genes of interest. Co-expression analysis based on Spearman's Correlation was performed using cBioPortal (http://www.cbioportal.org/) (25).

Survival analysis

Data from patients with HCC derived from TCGA were divided into two groups according to the expression of APOA2 (high or low, with an mRNA Z-score >0 or <0, respectively) and the expression of JUNB (high or low, with an mRNA Z-score >0 or <0, respectively). Survival curves were estimated by the Kaplan-Meier method and compared with the log-rank test.

Transcriptional regulation analysis

For each gene, interactions between proteins and their promoters and enhancers were obtained from r GeneHancer (26), a database of ChIP-seq data classify by inferred target genes. The interactions between TFs and their binding sites in the promoter and enhancer regions were supported by ChIP-sequencing. Subsequently, the TFs predicted to regulate both JUNB and APOA2 were used.

Results

Identification of significant differences in gene expression between in vivo and in vitro cells

Differential expression analysis was performed between the transcriptomes of in vivo tumor cells and two in vitro cell lines (Fig. S1). Top 2 cell population-specific marker gene expression is presented in Fig. S2A. Top 9 cell type-specific expressed genes are presented in Fig. S2B (HuH7 cells, Fig. S2Ba; P1T and P1C cells, Fig. S2Bb; P1B cells, Fig. S2Bc; and HuH1 cells; Fig. S2Bd). The specific expression was potentially due to the invasion of immune cells in the tumor and the immune response elicited by tumor cells (27). These significant differences suggested that the validation of the RNA-seq data should be a crucially important component for the in vitro analysis of the immune response in tumors. Cell population identifications of PBMCs are presented in Fig. S3A-D. Subsequently, an integrated comparison among tumor cells, cell lines and PBMCs was performed (Fig. 1). As expected, HuH1 and HuH7 cells clustered together in the analysis (Fig. 1A). Interestingly, the clustering analysis showed that a small number of epithelial cells were observed among the PBMCs (Fig. 1A-C). According to previous studies (28–30), these epithelial cells may be vascular endothelial cells, which have the potential to regulate the tumor cells. In addition to B cells and T cells, numerous types of mononuclear cells were detected among the PBMCs (Fig. 1C). A heat map was used to observe the expression of the differentially expressed genes (marker genes) in the comparisons (Fig. 1D). The expression of the marker genes in each group of single cells is presented in Fig. 1E.

Figure 1.

Single-cell transcriptomic data of PBMCs, HCC tumor cells and liver cancer cell lines integration. (A) Integration of PMBC and tumor samples. (B) Identification of cell populations by PCA. (C) Definition of cell types according to the marker genes. (D) Heatmap of the markers identified by PCA. (E) Comparison between expression of marker genes in PBMC and tumor samples. The color blue indicates the expression levels in PBMC samples, the color red in tumor samples. The size of each dot represents the percentage of the cells that expressed the corresponding gene. H1, HuH1; H7, HuH7; Mo, monocytes; NK, natural killer cells; P, epithelial cells; T, T cells; B, B cells; PBMC, peripheral blood mononuclear cells; PCA, principal component analysis; STIM, tumor samples; CTRL, PBMC samples; HCC, hepatocellular carcinoma.

Tumor-infiltrating immune cells exhibit a higher transcription level

In the analysis of differentially expressed genes, the gene expression levels between tumor cells and vascular epithelial cells in peripheral blood (Fig. 2A) were compared. In addition, the differentially expressed genes between T cells infiltrated in tumors and T cells in peripheral blood were examined (Fig. 2B). In the two comparisons, the number of the differentially expressed genes in T cells was significantly decreased compared with those in epithelial cells, while most of the differentially expressed genes in T cells were due to differences among cell types (Fig. 2C and D). Enrichment analysis of these differentially expressed genes showed that the differentially expressed genes in epithelial cells were significantly enriched in ‘inflammation’ and ‘immune response’ (Fig. 2E). The inhibition of some positive regulators of T cells, such as tumor necrosis factor (TNF), NF-κB inhibitor α (NFKBIA), Fos proto-oncogene (FOS), JUN and DEAD-box helicase 3 X-linked (DDX3X) may lead to the downregulation of the T cell receptor signaling pathway, and the inhibition of the hepatitis B and cellular oxidant detoxification pathways, which may be caused by the abnormal growth of tumor cells and the accumulation of stress-associated factors. In this condition, tumor cells may require more energy to sustain their proliferation and to adapt to the hypoxic micro-environment.

Figure 2.

Differentially expressed analysis in tumor cells and T cells. (A) A scatter plot of the gene expression correlation between tumor cells and vascular epithelial cells in the peripheral blood. (B) A scatter plot showing the gene expression correlation between T cells infiltrated in the tumor and T cells in the peripheral blood. (C) A heatmap of differentially expressed genes in the comparisons. (D) A Venn diagram of the number of the differential genes with |log2Fold-Change|>2 in the comparisons. (E) Enriched Gene Ontology terms and Kyoto Encyclopedia of Genes and Genomes pathways of these differential genes. UMI, unique molecular identifier; PBMC, peripheral blood mononuclear cells.

JUNB may serve a crucial role in HCC tumor immune microenvironment

JUNB, which was reported as highly expressed in cancer cells in previous studies (22,31–33), was observed to be downregulated in numerous tumor cells, although it was observed to be highly expressed in certain tumor cell lines (Fig. 3A). Additionally, all PBMCs presented high expression levels of JUNB (Fig. 3B). Since JUNB is associated with lipid metabolism, high expression of JUNB was expected in activated cells. S100A4, which is associated with the development and progression of the tumor and immune infiltration (34,35), exhibited high expression in in vivo epithelial cells and T cells, but low expression in in vitro cell lines (Fig. 3D and E). Although the indirect interaction between JUNB and S100A4 was found to function through annexin A2 (ANXA2; Fig. S3E), a positive association between ANXA2, JUNB and S100A4 was detected in patients with HCC from the TCGA dataset (Fig. S3F). In addition, the functional roles of ANXA2, JUNB and S100A4 in tumor cells were found to be associated with T cells (14,36,37), suggesting that the function of JUNB in HCC tumor cells may be associated with the interaction between these genes in the tumor immune microenvironment. By combining the data from PBMCs and tumor samples, JUNB and S1004A were found to be decreased in epithelial cells (Fig. 3C and F). However, the expression level of S100A4 was found to be increased in tumor-associated T cells, while the expression of JUNB was decreased. JUNB was identified to be downregulated in various tumors in the previous studies (31).

Figure 3.

Violin plots of the gene expression levels of JUNB and S100A4 in different cell populations. (A) Gene expression of JUNB in tumor samples. (B) Gene expression of JUNB in PBMC samples. (C) Gene expression of JUNB in integrated samples. (D) Gene expression of S1004A in tumor samples. (E) Gene expression of S1004A in PBMC samples. (F) Gene expression of S1004A in integrated samples. PBMC, peripheral blood mononuclear cells; UMI, unique molecular identifier.

Validation of the key role of JUNB in an independent dataset

The potential role of JUNB was further investigated by integrating a cohort of 360 HCC patients from the TCGA database. APOA2 was found to be significantly negatively associated with JUNB in the dataset (P=1.094×10−5; Fig. 4A). Among the 360 HCC samples, >20% of the samples showed APOA2 mRNA upregulation, which may be explained by the low expression of JUNB in HCC samples (Fig. 4B). The high expression of JUNB was previously predicted as a poor survival indicator in patients with tumors (32,33), while APOA2 was found to serve a role in the development and progression of the tumor through the peroxisome proliferator activated receptor α (PPARα) pathway (38). Moreover, >50% of the TFs that interact with the promoters and enhancers of these two genes were found to regulate both APOA2 and JUNB (Fig. 4C). Further data enrichment analysis showed that the pathways associated with the TFs regulating APOA2 and JUNB were involved in the development and progression of HCC and immune response, which may represent a potential novel mechanism underlying HCC (Fig. 4D), although further experimental evidence is required to test this hypothesis. Although survival analysis showed a slightly different result when analyzing the overall survival time in two groups of patients (P=0.11), a significant difference was observed when analyzing only longer survival time (>15 months; Fig. 4E), indicating that JUNB and APOA2 may play a key role in improving the survival time in patients with HCC.

Figure 4.

Validation of JUNB expression in the TCGA dataset. (A) A scatter plot of the gene expression correlation of APOA2 and JUNB. (B) Mutation and expression level of APOA2 and JUNB in the HCC dataset from TCGA. (C) A Venn diagram of the number of the regulatory elements of APOA2 and JUNB detected in GeneHancer. (D) TFs that coregulated APOA2 and JUNB. Color of the edge indicates the confidence level of the interaction between two proteins. The nodes, in red, indicate the PPAR signaling pathway associated proteins. (E) A survival curve for patients from TCGA with high/low expression of APOA2 and low/high expression of JUNB. TF, transcription factor; TCGA, The Cancer Genome Atlas; APOA2, apolipoprotein A2; PPAR, peroxisome proliferator activated receptor; HCC, hepatocellular carcinoma; JUNB, JunB proto-oncogene.

Discussion

Increased understanding of tumor-host interactions has accelerated the development of novel cancer immunotherapies. In addition, drug resistance in biomarker therapy and immunotherapy have recently been investigated (21). Despite their success, immune-checkpoint inhibitors present certain limitations. Therefore, in addition to improving the treatment of patients presenting with tumors at an advanced stage, the identification of driver biomarkers may contribute to immunotherapy in an early clinical stage (39). Since immune-checkpoint inhibitors could be used in HCC treatment, combining molecular targeted therapy with immunotherapy has become a therapeutic method to stimulate the immune response. The previously described mechanisms underlying tumor development based on cell lines studies have been found to be unreliable, particularly for immune response and immunotherapy-related studies. In order to investigate the role of the immune response in tumor progression, the tumor microenvironment and the balance between tumor cells and immune cells must be considered.

The present results suggested that the inhibition of JUNB may be a key indicator of the regulation of the APOA2-associated PPARα pathway in HCC (31). APOA2 is a well-known member of the apolipoprotein family (40), which is functionally involved in triglyceride, fatty acid and glucose metabolism. This gene family has been previously reported to be overexpressed in HCC and regulated by the PPARα pathway (41). In addition, together with the co-expression of APOA2 and JUNB observed in the TCGA dataset and the transcriptional regulation analysis, a potential regulation of APOA2 and JUNB by the PPARα pathway was identified. Finally, the present survival analysis in HCC patients suggested that the investigation of JUNB may facilitate the identification of novel therapeutic targets for HCC patients. The inhibition of some positive regulators of T cells, such as TNF, NFKBIA, FOS, JUN and DDX3X may lead to the downregulation of the T cell receptor signaling pathway (42–44), affecting the hepatitis B and cellular oxidant detoxification response, which may be caused by the abnormal growth of tumor cells and the accumulation of stress-associated factors. In this condition, tumor cells may require more energy to sustain their proliferation and to adapt to the hypoxic micro-environment.

Although the present results could potentially contribute to the development of novel therapeutics to treat patients with HCC, the lack of experimental validations on mRNA levels are still the main limitations for the study. In addition, investigating the expression levels of the products of candidate mRNAs, the proteins, is equally critical to validate the results, as they play a central role in biological processes. Furthermore, validations using both in vivo and in vitro models should be further investigated, such as ChIP to detect the potential regulatory regions on JUNB, Junb defective mice to study the important functions of Junb.