Comprehensive characterization of enhancer RNA in hepatocellular carcinoma reveals three immune subtypes with implications for immunotherapy.

Hepatocellular carcinoma (HCC) is highly heterogeneous. Molecular subtyping for guiding immunotherapy is warranted. Previous studies have indicated that enhancer RNAs (eRNAs) are involved in tumor heterogeneity and immune infiltration. However, the eRNA landscape and its correlation with immune infiltration in HCC remain unknown. Here we first revealed the genome-wide eRNA landscape in two HCC cohorts. Then we divided individuals with HCC into three immune-related clusters (C1, C2, and C3) based on eRNA expression profiles. The prognosis, biological properties, immune infiltration, clinical features, genomic features, and drug response were analyzed. C1 was enriched in immune infiltration and potentially sensitive to immune checkpoint inhibitors (ICIs). C2 displayed features of immune depletion, high proliferation activity, malignant clinical features, and the worst prognosis. C2 may benefit from targeted therapy. C3 presented moderate immune infiltration, metabolism-related signatures, and the best prognosis. Transarterial chemoembolization (TACE) may be effective for C3. Finally, we constructed a 51-eRNA classifier for subtype prediction and validated its efficacy in The Cancer Genome Atlas (TCGA) cohort and Sun Yat-sen University Cancer Center (SYSUCC) cohort. Our results provide a novel method for immune classification of HCC, shed new light on tumor heterogeneity, and may aid in HCC immunotherapy.

Introduction

Liver cancer is the second leading cause of cancer-related death worldwide. Hepatocellular carcinoma (HCC) is the most prevalent pathological type of liver cancer. Although several recent improvements have been made in comprehensive therapy of HCC, the 5-year survival of individuals with HCC still remains poor. HCC is a highly heterogeneous disease. Molecular subtyping has important implications when devising individualized therapeutic strategies for individuals with HCC.

Enhancers are essential regulatory elements of DNA that activate target gene transcription., Tens of thousands of enhancers have been identified in the Encyclopedia of DNA Elements (ENCODE). Although enhancers generally regulate gene expression in cis, they also exert long-range effects, which are independent of position and orientation when activated by chromosomal rearrangement or specific transcription factors.8, 9, 10 Recently, with the development of RNA sequencing technology, it has been reported that enhancers also function as noncoding RNA transcription units to generate enhancer RNAs (eRNAs), which can be used to infer the activities of enhancer loci. As an important component of the transcriptome and epigenome, accumulating evidence shows that eRNAs are widely involved in human cancer by modulating chromatin structure or interacting with transcriptional regulators.12, 13, 14 Qin et al. performed molecular subtyping of lung adenocarcinoma based on functional eRNAs, which identified distinct immune signatures among different subtypes and might be instructive for cancer immunotherapy. However, the genome-wide eRNA landscape and its correlation with the immune microenvironment in HCC remain largely undefined.

In this study, we illustrated the genome-wide landscape of eRNAs in HCC by analyzing total RNA sequencing data from the GEO database. We established a novel immune-related HCC subtyping strategy based on the expression profile of eRNAs, and the prognosis, immune infiltration, biological properties, clinical features, genomic features and drug responses of different subtypes were analyzed. Finally, we generated a 51-eRNA classifier for subtype prediction. The efficacy of the classifier was validated in HCC cohorts from The Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) and Sun Yat-sen University Cancer Center (SYSUCC).

Results

Landscape of eRNAs in HCC

The main research processes of this study are summarized in Figure 1. We first explored the genomic distribution of eRNAs and found that they were relatively uniform across all chromosomes except chromosome 13 (Figure 2A). The length of most eRNAs was around 300 bp in two cohorts (Figure 2B). The average expression level of global eRNA was significantly upregulated in HCC compared with that in normal liver (Figure 2C), implying that eRNAs may be involved in HCC carcinogenesis. Next we analyzed the average expression level of global eRNA in different groups based on age, sex, α-fetoprotein (AFP) level, tumor size, Tumor Node Metastasis (TNM) stage and hepatitis viral infection status. No significant difference in the average expression level of global eRNA was observed between different age, sex, AFP level, and tumor size groups (Figures S1A–S1D), except for remarkably downregulated expression of global eRNAs in hepatitis C virus (HCV)-infected samples versus hepatitis B virus (HBV)-infected samples in the Yoon et al. cohort (Figure S1F), which suggested a distinct regulation pattern because of different virus infections. eRNAs tended to be highly expressed in advanced HCC tissues in the Candia et al. cohort, although it was not statistically significant (p = 0.28; Figure S1E).

Figure 1

Overview of the study workflow

Figure 2

Landscape of eRNA in HCC

(A) The genomic distribution of eRNA in the Candia et al. (left) and Yoon et al. cohorts (right). Pink dots, locations; purple peaks, density of eRNA distribution. (B) Length distribution of eRNAs in two cohorts. (C) Average expression of global eRNAs in HCC and normal liver tissues in two cohorts. (D) Volcano plot of DEeRNAs in HCC and normal liver tissue in two cohorts. Red, upregulated; white, not significantly different; blue, downregulated. (E) The distribution of DEeRNAs on 23 chromosomes in two cohorts. Orange dots, location of upregulated eRNAs; blue dots, location of downregulated eRNAs; peaks, density of eRNA distribution. (F) Venn diagram of DEeRNAs in two cohorts. (G) Functional enrichment analysis of upregulated eRNAs in two cohorts. Boxplot “boxes” indicate the first, second, and third quartiles of the data.∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, ∗∗∗∗p < 0.0001.

Differentially expressed eRNAs (DEeRNAs) in HCC and normal liver were identified (Figure 2D). There were 68 upregulated eRNAs and 55 downregulated eRNAs in both cohorts (Figure 2F). The genomic distribution of DEeRNAs is presented in Figure 2E. Functional enrichment analyses of co-expressed protein-coding genes (PCGs) of DEeRNAs were performed. We focused on eRNAs upregulated in HCC and found that PCGs correlated to these eRNAs were mainly involved in cell cycle, DNA replication, mismatch repair, and several metabolism-related pathways, indicating that these eRNAs may play essential roles in HCC carcinogenesis (Figure 2G; Table S2).

Identification of three HCC subtypes based on expression profiles of eRNAs

To explore HCC molecular subtyping based on expression profiles of eRNAs, we conducted non-negative matrix factorization (NMF) analyses in two training cohorts (the Candia et al. and Yoon et al. cohorts). The cophenetic correlation coefficients were calculated to decide the appropriate number of subtypes, and k = 3 was chosen (Figures 3A and 3D). We classified individuals with HCC into three subtypes: clusters 1, 2, and 3 (C1, C2, and C3, respectively) (Figures 3B and 3E). Principal-component analysis (PCA) showed that individuals in the three subtypes were well separated (Figures 3C and 3F). SubMap analyses confirmed that C1, C2, and C3 in the Candia et al. cohort were highly associated with corresponding subtypes in the Yoon et al. cohort (Figure 3G). Significant differences between subtypes were observed in overall survival (OS) or disease-free survival (DFS). C3 had a survival advantage over C1 and C2 (OS in the Candia et al. cohort, C3 versus C1, p = 0.029, C3 versus C2, p < 0.001; DFS in the Yoon et al. cohort, C3 versus C1, p = 0.025, C3 versus C2, p < 0.001) (Figures 3H and 3I). The subtype-specific eRNAs are shown in Figures S2A and S2B and Table S4.

Figure 3

Identification of three HCC subtypes based on expression profiles of eRNAs

(A) Cophenetic correlation coefficient for k = 2–7 for the Candia et al. cohort. (B) NMF subtyping based on eRNA expression profile in the Candia et al. cohort. (C) PCA of eRNAs to divide individuals with HCC into three subtypes in the Candia et al. cohort. (D) Cophenetic correlation coefficient for k = 2–7 for the Yoon et al. cohort. (E) NMF subtyping based on eRNA expression profile in the Yoon et al. cohort. (F) PCA of eRNAs to divide individuals with HCC into three subtypes in the Yoon et al. cohort. (G) SubMap analysis of GEPs between the HCC subtypes in the Candia et al. and Yoon et al. cohorts. (H) OS of three subtypes in the Candia et al. cohort. (I) DFS of three subtypes in the Yoon et al. cohort.

Prognostic significance of the eRNA subtyping system in HCC

We first evaluated the prognostic significance of the eRNA subtyping system in the Candia et al. and Yoon et al. cohorts. Univariate and multivariate Cox regression identified C3 as an independent prognostic factor in the Candia et al. and Yoon et al. cohorts (Figures 4A and 4D). Then we evaluated the efficacy of the eRNA subtyping system in prognosis prediction through time-dependent receiver operating characteristic (ROC) analyses. The area under the ROC curve (AUC) for 1-, 3- and 5-year OS was 0.8, 0.71, and 0.64 in the Candia et al. cohort and 0.71, 0.72, and 0.79 for 1-, 2-, and 3-year DFS in the Yoon et al. cohort (Figures 4B and 4E). We also compared the predictive ability of eRNA subtype with other clinical parameters. The results showed that the predictive ability of the eRNA subtyping system was superior to the predictive ability of TNM staging system, serum AFP level, tumor size, or cirrhosis (Figures 4C and 4F).

Figure 4

Prognostic significance of the eRNA subtyping system in HCC

(A) Univariate and multivariate Cox regression analyses of OS in the Candia et al. cohort. (B) Time-dependent ROC analysis of eRNA subtyping in the Candia et al. cohorts. (C) The sensitivity and specificity of eRNA subtyping, TNM stage, and cirrhosis in prognostic prediction were estimated by time-dependent ROC curves in the Candia et al. cohort. (D) Univariate and multivariate Cox regression analyses of DFS in the Yoon et al. cohort. (E) Time-dependent ROC analysis of eRNA subtyping in the Yoon et al. cohort. (F) The sensitivity and specificity of eRNA subtyping, TNM stage, tumor size, and AFP level in prognostic prediction were estimated by time-dependent ROC curves in the Yoon et al. cohort.

Biological properties of the three HCC subtypes

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and cancer hallmark gene set enrichment analyses were utilized to distinguish the characteristic pathways related to the three subtypes. As shown in Figures 5A and 5B, pathways involved in metastasis, immune infiltration, and tumorigenesis were significantly enriched in C1, indicating that C1 was correlated to immune infiltration. C2 was inactive in immune-related pathways but had a remarkable enrichment in proliferation-related pathways, including G2M checkpoint, DNA replication, mismatch repair, E2F targets, and cell cycle. C3 presented moderate activity in immune-related pathways and was activated in metabolism-related pathways, including adipogenesis, bile acid metabolism, and drug metabolism, suggesting that C3 was a well-differentiated subtype of HCC. Generally, C1 had a more activated immune microenvironment compared with C2 and C3.

Figure 5

Functional enrichment analyses of three HCC subtypes

(A and B) Heatmaps of biological properties enriched in three subtypes in the Candia et al. (A) and Yoon et al. (B) cohorts, analyzed by ssGSEA. High and low ssGESA scores are represented in red and blue, respectively. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, ∗∗∗∗p < 0.0001.

Characterization of immune infiltration in the three HCC subtypes

Pathway enrichment analyses revealed that the extent of immune infiltration was significantly different in the three subtypes. We analyzed the immune and stromal scores according to the estimation of stromal and immune cells in malignant tumours using expression data (ESTIMATE) algorithm (Figures 6A and 6B). C1 displayed markedly higher immune and stromal scores than C2 and C3 in the Candia et al. cohort. A similar trend was observed in the Yoon et al. cohort, which was in line with the results of functional enrichment analyses. C1 also exhibited higher cytolytic activity scores than C2 and C3, indicating activation of the tumor microenvironment (TME) in C1 (Figure 6C).

Figure 6

Correlation of the HCC subtypes with immune infiltration

(A and B) The immune scores and stromal scores of HCC subtypes in the Candia et al. and Yoon et al. cohorts. (C) The cytolytic activity scores of HCC subtypes in the Candia et al. and Yoon et al. cohorts. (D) Expression level of 15 immune checkpoint genes in three HCC subtypes in two cohorts. (E) Enrichment of immune and stromal cells in three HCC subtypes in two cohorts, analyzed by TIMER, MCP-counter, and xCell, respectively. Boxplot “boxes” indicate the first, second, and third quartiles of the data. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, ∗∗∗∗p < 0.0001.

We also evaluated the expression level of 15 immune checkpoint genes that may become potential therapeutic targets for HCC in the three subtypes. C1 showed the highest expression of 13 genes (except for LAG3 and OX40L) among the three subtypes, suggesting that individuals in C1 may benefit from immune checkpoint inhibitors (ICIs). Both cohorts exhibited similar results (Figure 6D).

Then we analyzed immune infiltration cells in the three HCC subtypes through tumor immune estimation resource (TIMER), microenvironment cell populations (MCP)-counter, and xCell. C1 showed significant enrichment of immune cells, including B cells, CD8+ T cells, and dendritic cells. C1 and C3 displayed higher enrichment of endothelial cells and cancer-associated fibroblast cells. Consistent with the immune and stromal scores, C2 showed the lowest activity in immune cells and stromal cells (Figure 6E). Our eRNA subtyping system distinguished three immune-related HCC subtypes. C1 was defined as an immune-enriched subtype and C2 as an immune-depleted subtype. C3 was moderately infiltrated by immune cells.

Clinical characteristics of the three HCC subtypes

We next investigated clinical features of the three subtypes in the Candia et al. and Yoon et al. cohorts. In the Candia et al. cohort, C2 had a higher level of AFP (Figure 7A). Similarly, C2 also had more aggressive clinical features in the Yoon et al. cohort, such as larger tumor size and increased relapse rate (Figure 7B). These results confirmed that C2 was correlated with poor prognosis. Other clinical characteristics, such as age, sex, virus infection, tumor numbers, and Child-Pugh score, showed no significant difference among the three subtypes.

Figure 7

Clinical characteristics of HCC subtypes in two cohorts

(A and B) Clinical characteristics of HCC subtypes in the Candia et al. and Yoon et al. cohorts. (C and D) Association of eRNA-based subtypes with previous HCC classifications in the Candia et al. and Yoon et al. cohorts. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, ∗∗∗∗p < 0.0001.

Considering that several HCC classifications have been reported so far, we estimated the correlation between eRNA subtype and previous classifications, including the Boyault et al. subtype, Hoshida et al. subtype, Lee et al. subtype, and Roessler et al. subtype (Figures 7C and 7D). In the Candia et al. cohort, the C1 subtype was significantly linked to the Boyault et al. G4 class, the Hoshida et al. S1 class, the Lee et al. A class, and the Roessler et al. high-risk group (p < 0.001). The C2 subtype was strongly correlated with the Boyault et al. G2/3 class, the Hoshida et al. S2 class, the Lee et al. A class, and the Roessler et al. high-risk group (p < 0.001). The C3 subtype was associated with the Boyault et al. G4/6 class, the Hoshida et al. S3 class, the Lee et al. B class, and the Roessler et al. low-risk group (p < 0.001). Similarly, in the Yoon et al. cohort, the C1 subtype was significantly correlated with the Boyault et al. G4 class, the Hoshida et al. S1 class, the Lee et al. A class, and the Roessler et al. high-risk group (p < 0.001). The C2 subtype was strongly linked to the Boyault et al. G4/5 class, the Hoshida et al. S2 class, and the Roessler et al. high-risk group (p < 0.001). The C3 subtype was associated with the Boyault et al. G4 class, the Hoshida et al. S3 class, the Lee et al. B class, and the Roessler et al. low-risk group (p < 0.001). These results revealed that eRNA subtype had great consistency with previous studies.

Genomic features of different immune subtypes

To investigate whether gene somatic mutation was distinct among the three subtypes, we focused on genes involved in cell cycle regulation, Wnt/β-catenin signaling, hepatic differentiation, and chromatin modification in the Candia et al. cohort (Figures 8A and 8B). Generally, the C1 subtype had the lowest gene mutation burden among the three subgroups (Figure 8C). The mutation frequency of cell-cycle-related genes was higher in C2 (53%) than C1 (22%) and C3 (46%), which was consisted with the proliferative feature of C2. Although the mutation frequency of cell-cycle-related genes was also common in C3, the function enrichment result and the correlation analysis with other HCC classifications did not support C3 as a high-proliferation subtype. Albumin (ALB) and apolipoprotein B (APOB) mutation were most common in the C3 subtype (50%) compared with C1 (5.5%) and C2 (21%), in agreement with the dysregulation of cellular metabolism in C3. CTNNB1 mutation was most frequent in the C3 subtype (39.3%), suggesting that the Wnt/β-catenin pathway may be dysregulated in C3. No significant difference in chromatin-modification-related gene mutations was observed among the three subtypes (Figure 8A).

Figure 8

Genomic features of three HCC subtypes

(A) Gene mutation landscape of HCC subtypes in the Candia et al. cohort. (B) Mutation frequency of genes in the Candia et al. cohort. (C) Total TMB of HCC subtypes in the Candia et al. cohort. Boxplot “boxes” indicate the first, second, and third quartiles of the data.∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, ∗∗∗∗p < 0.0001.

Potential response to targeted therapy, immunotherapy, and transarterial chemoembolization (TACE) in the three HCC subtypes

Considering that the immune microenvironment is closely related to drug response in cancer treatment, we explored whether eRNA-based classification could be instructive for developing an individualized therapeutic strategy in clinical practice. We compared the gene expression profiles (GEPs) of the three immune subtypes with HCC cell lines that had drug sensitivity data. For targeted agents (sorafenib and cabozantinib), C2 showed a significant correlation with the sorafenib and cabozantinib response group (sorafenib, p = 0.003; cabozantinib, p = 0.015) in the Yoon et al. cohort (Figure 9A). The same tendency was observed in the Candia et al. cohort, but the results were not all statistically significant (sorafenib, p = 0.015; cabozantinib, p = 0.053) (Figure 9D). On the contrary, C1 was strongly associated with the sorafenib- and cabozantinib-resistant group (Candia et al. cohort, sorafenib, p = 0.004, cabozantinib, p = 0.018; Yoon et al. cohort, sorafenib, p = 0.002, cabozantinib, p = 0.011) (Figures 9A and 9D). These results showed that C2 was more likely to benefit from sorafenib and cabozantinib treatment.

Figure 9

Correlation of the HCC subtypes with the response to targeted therapy, immunotherapy, and TACE

(A–C) SubMap correlation analysis between HCC subtypes in the Candia et al. cohort and samples with different sensitivities to targeted therapy (A), ICIs (B), and TACE (C). (D–F) SubMap correlation analysis between HCC subtypes in the Yoon et al. cohort and samples with different sensitivities to targeted therapy (D), ICIs (E), and TACE (F).

C1 was more likely to respond to the programmed cell death protein 1 (PD-1) inhibitor in both cohorts (Candia et al. cohort, p = 0.014; Yoon et al. cohort, p = 0.033) (Figures 9B and 9E). The same tendency was observed for cytotoxic T lymphocyte-associated protein 4 (CTLA-4) inhibitor response, but the results were not all statistically significant (Candia et al. cohort, p = 0.131; Yoon et al. cohort, p = 0.033). However, C2 seemed to be immunotherapy resistant because of immune depletion features in the Candia et al. cohort (PD-1 inhibitor, p = 0.008; CTLA-4 inhibitor, p = 0.045) (Figure 9B).

TACE is a widely used treatment for HCC., We analyzed the potential TACE response status in the three HCC subtypes and found that the GEP of C3 was significantly correlated with the TACE response group (Candia et al. cohort, p = 0.001; Yoon et al. cohort, p = 0.004) (Figures 9C and 9F), indicating that C3 tends to benefit from TACE treatment.

Construction and validation of a 51-eRNA classifier for eRNA subtype identification

For subtype prediction, we performed differential expression analysis to construct a subtype-specific eRNA classifier. After integrating subtype-specific eRNAs in the Candia et al. and Yoon et al. cohorts, a total of 51 eRNAs, called the 51-eRNA classifier, were screened out (Figure S2C; Table S5). High accordance was observed between subtypes predicted by the 51-eRNA classifier based on NTP and original classification based on NMF, with concordance of 80.77% in C1, 72.72% in C2, and 95.45% in C3 in the Candia et al. cohort and 88.24% in C1, 100% in C2, and 73.68% in C3 in the Yoon et al. cohort (Figure 10A).

Figure 10

Construction and validation of the subtype-specific eRNA classifier

(A) Concordance of HCC subtypes between the prediction based on the 51-eRNA classifier and original classification based on NMF in two cohorts. (B) OS of three HCC subtypes predicted by the 51-eRNA classifier in the TCGA-LIHC cohort.

(C and D) The immune scores, stromal scores, and cytolytic activity scores of HCC subtypes predicted by the 51-eRNA classifier in the TCGA-LIHC cohort (C) and SYSUCC cohort (D). (E and F) Enrichment of immune and stromal cells in three HCC subtypes predicted by the 51-eRNA classifier in the TCGA-LIHC cohort (E) and SYSUCC cohort (F), analyzed by TIMER, MCP, and Xcell, respectively. (G) Representative images of IHC-stained slides with CD8 and ki67 antibodies in three subtypes in the SYSUCC cohort. Scale bars, 100 μm. Boxplot “boxes” indicate the first, second, and third quartiles of the data. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, ∗∗∗∗p < 0.0001.

To test the efficacy of the subtype-specific eRNA classifier, we divided the samples in the two validation cohorts into three subtypes based on the 51-eRNA classifier using the NTP algorithm. The prognosis, immune infiltration, biological properties, genomic features, and drug responses among the three subgroups were analyzed. As expected, C1-predicted displayed the highest levels of immune cell infiltration and most active immune-related pathways. C2-predicted was still characterized by least immune cell infiltration, high proliferation, and poor prognosis. C2-predicted showed a higher copy number alteration (CNA) burden than the C1- and C3-predicted subtypes in the TCGA-LIHC cohort. C3-predicted presented moderate immune cell infiltration and metabolism signatures (Figures 10B–10F, S3, and S4). Similar results were observed in the SYSUCC cohort. Immunohistochemistry (IHC) staining demonstrated the highest CD8+ T cell infiltration in the C1 subtype (Figure 10G). C2 displayed the most abundant ki67-positive cells, consistent with high proliferation in C2 (Figure 10G). These results showed that our eRNA subtypes were highly reproducible and that the 51-eRNA classifier was reliable.

In accordance with the result of the Kaplan-Meier analysis (Figure 10B), univariate Cox regression revealed that C3 had a lower survival risk than C2 in the TCGA-LIHC cohort (Figure S5A). The AUC of the 51-eRNA classifier for 1-, 3-, and 5-year OS was 0.64, 0.60, and 0.57 in the TCGA-LIHC cohort (Figure S5B). The efficacy of the eRNA subtyping system was shown to be superior to the TNM staging system and vascular invasion in prognosis evaluation (Figure S5C).

The efficacy of the classifier in prediction treatment responses was also validated in the TCGA-LIHC and SYSUCC cohorts. C1-predicted was sensitive to ICIs (TCGA-LIHC cohort, PD1 inhibitor, p = 0.003, CTLA4 inhibitor, p = 0.049 [Figure 11B]; SYSUCC cohort, PD1 inhibitor, p = 0.001 [Figure 11E]). C2-predicted still showed a significant correlation with the sorafenib- and cabozantinib-response group (TCGA-LIHC cohort, sorafenib, p = 0.001, cabozantinib, p = 0.001 [Figure 11A]; SYSUCC cohort, sorafenib, p = 0.003, cabozantinib, p = 0.040 [Figure 11D]). C3-predicted was likely to benefit from TACE (TCGA-LIHC cohort, p = 0.001 [Figure 11C]; SYSUCC cohort, p = 0.002 [Figure 11F]). The 51-eRNA classifier is highly valuable for predicting therapeutic responses of individuals with HCC.

Figure 11

Predictive value of the 51-eRNA classifier in response to targeted therapy, immunotherapy, and TACE

(A–C) SubMap correlation analysis between HCC subtypes in the TCGA-LIHC cohort and samples with different sensitivities to targeted therapy (A), ICIs (B), and TACE (C). (D–F) SubMap correlation analysis between HCC subtypes in the SYSUCC cohort and samples with different sensitivities to targeted therapy (D), ICIs (E), and TACE (F).

A summary of the three immune-related HCC subtypes is shown in Figure 12.

Figure 12

Overview of three HCC subtypes based on eRNA expression profile

Discussion

HCC is a highly heterogeneous disease. Several staging systems for HCC (e.g., TNM and the Barcelona clinic liver cancer staging system) are generally based on clinical information and lack a systemic overview of molecular features, which is insufficient for guiding individualized therapy, such as immunotherapy and target therapy., Most molecular models of HCC are based on gene expression and mutation.,,24, 25, 26 Few studies have been performed to explore epigenetic signatures during tumorigenesis.

eRNAs are important parts of the epigenetic regulation system transcribed by activated enhancers., Accumulated evidence indicates that eRNAs are strongly associated with cancer development. In this study, we characterized the eRNA landscape of HCC and uncovered its important role in the cell cycle, extracellular matrix (ECM)-receptor interaction, DNA replication, mismatch repair, and metabolism. We constructed a novel three-subtype molecular classification system for HCC based on eRNA expression profiles. The prognosis, immune signature, clinical features, biological properties, and drug responses were distinct among the three subtypes. C1 was distinguished by enriched immune infiltration and high expression of immune-checkpoint-related genes, which implied that C1 samples were “hot” tumors and sensitive to ICIs. C2 exhibited a lack of immune infiltration, high proliferation activity, and poor prognosis. C2 might benefit from targeted therapy. C3 was defined as a moderate immune subtype with metabolism-related signatures, low AFP level, and good prognosis. TACE was predicted to be effective for C3.

The tumor immune microenvironment has been recognized as a key factor in HCC biological behavior and drug response. Immune classification has the potential to predict prognosis and treatment responses for individuals with HCC. In recent years, immunotherapy has achieved inspiring outcomes for individuals with HCC.32, 33, 34 In this study, we found that C1 represented a small group of individuals with HCC who might be suitable for immunotherapy. The high immune score and enriched immune infiltration of C1 are in line with the notion that tumors with abundant immune cell infiltration are more likely to respond to immunotherapy. C1 accounts for about 30% of individuals in every HCC cohort, which is in agreement with the objective response rate of ICIs in clinical trials.,, C2 presented deficiency in immune infiltration cells, suggesting that immunotherapy may be ineffective for C2. The enrichment of TP53 mutations and high proliferation of C2 may contribute to its resistance to ICIs. Previous studies have revealed that sorafenib responders were more likely to be highly proliferative and ICI resistant, which is consist with the C2 subtype in this study.

Previous research has identified several HCC subgroups associated with recurrence or OS.,,24, 25, 26 In our study, C1 showed great enrichment of metastasis-related pathways, which is associated with the Roessler et al. high-risk subgroup. Notch and transforming growth factor β (TGF-β) signaling were aberrantly activated in C1, which indicates a strong correlation with the Hoshida et al. S1 subgroup. C2 was distinguished by high proliferation and poor prognosis, which is similar to Boyault et al. subtype G2. Individuals in C2 had the highest frequency of TP53 mutation. The TP53 mutation is strongly associated with cell proliferation and chromosome instability, which may partially result in the malignant features of C2. The high AFP level of C2 is in concordance with the “progenitor”-like feature of the Hoshida et al. S2 subgroup. However, a controversial result was observed in the correlation analyses between the Lee et al. subtype and eRNA subtype in different HCC cohorts, which may result from different genetic background and etiology. C3 is strongly linked to the Hoshida et al. S3 subgroup. The good prognosis of C3 is also concordant with Lee et al. subtype B and the Roessler et al. low-risk subgroup. eRNA-based subtyping showed high consistency with published HCC subgroups, which implied that eRNAs were involved in HCC heterogeneity. Dissecting tumor heterogeneity from the perspective of eRNA expression is effective.

Development of HCC is in part due to multiple complex mutation processes. Aberrant Wnt/β-catenin signaling pathway activation because of CTNNB1 mutation has been observed in more than 20% of individuals with HCC., The highest mutant frequency of CTNNB1 was observed in the C3 subgroup (39%) compared with C1 (11.1%) and C2 (10%), suggesting dysregulation of the Wnt/β-catenin pathway in C3. It has been reported that activation of β-catenin signaling can lead to reduced immune infiltration and resistance to ICIs, which is consistent with our results showing that individuals in C3 may not benefit from immunotherapy. However, Wnt signaling pathway-targeted inhibitors may be effective for individuals in C3. A large proportion of C3 individuals carried ALB and APOB mutations, implying a highly differentiated status of C3, which may partially account for the good prognosis. Generally, a low tumor mutation burden (TMB) is associated with a poor response to ICIs.44, 45, 46 However, although the TMB of C1 was lower than those of the other two subtypes, C1 tended to respond to ICIs. C2 and C3 turned out to be insensitive to ICIs, partially because chronic antigen exposure because of high TMB may lead to dysfunction and exhaustion of T cells in C2 and C3. There is no significant difference in neoantigens among three subtypes in the TCGA cohort, suggesting that TMB is not the best predictor of immune infiltration and immunotherapy efficacy in HCC. Reduced immune-mediated cytotoxic and proinflammatory activity in high broad CNA-burden tumors had been observed before and was further proven in the C2 subtype. This finding was likely due to alterations in antigen-presenting machinery and decreased rates of neoantigenic mutations. High CNA burden was also linked to high proliferation, TP53 dysfunction, and DNA repair.

Targeted therapy and immune checkpoint blockade therapy epitomize new precision treatments for HCC. Sorafenib is a multitargeted kinase inhibitor and currently recommended as the first-line therapy of advanced HCC.50, 51, 52, 53 Another inhibitor of multiple receptor tyrosine kinases, cabozantinib, is a treatment option for individuals with HCC previously treated with sorafenib. As for immunotherapy, PD-1 inhibitors and CTLA-4 inhibitors also have promising applications in HCC.,, However, the low response rate and lack of efficacy predictor of ICIs still remain huge challenges. Here we integrated subtype-specific eRNAs into the Candia et al. and Yoon et al. cohorts to obtain a 51-eRNA classifier with the aim of classifying HCC subgroups and predicting responses to different therapies. The validation result showed that the C1, C2, and C3 subtypes were sensitive to ICIs, targeted therapy, and TACE, respectively. The high concordance of classifiers among different cohorts indicated that the eRNA classifier could reproducibly determine HCC classification and guide individualized therapy.

Our study has the following strengths. First, our study provides a new approach for HCC classification based on the expression profile of eRNAs. We utilized total RNA sequencing data to reveal the eRNA landscape in HCC, which is more credible than poly(A) enrichment sequencing. Most importantly, we established and validated a 51-eRNA classifier to predicting prognosis and treatment responses. Nevertheless, the study has some limitations. First, the sample size was relatively small because of the limited numbers of available total RNA sequencing cohorts. Second, comprehensive analyses together with other omics (e.g., proteomics and metabolomics) are needed in the future. Third, the biological roles of eRNAs in HCC carcinogenesis need further exploration.

Our results provide a novel method for immune classification of HCC, shed new light on tumor heterogeneity, and may aid in HCC immunotherapy.

Materials and methods

RNA sequencing of HCC samples

The institutional review board of SYSUCC approved this study, and written informed consent was obtained from all individuals. Baseline information of individuals in the SYSUCC cohort is shown in Table S1. A total of 14 fresh HCC tumor tissue samples were collected immediately after liver biopsy and stored in RNA stabilization solution (Invitrogen, USA) at −80°C. For preparing the transcriptome library, total RNA was extracted and treated with a rRNA removal kit (Invitrogen). After RNA integrity detection, cDNA synthesis, end repair, A-base addition, and ligation of index adapters were performed. Paired-end libraries were sequenced using the HiSeq X Ten platform.

Data collection and RNA sequencing analysis

Two groups of total RNA sequencing (RNA-seq) data from the Yoon et al. (GSE148355) and Candia et al. cohorts (GSE144269) were included in our study. Baseline information of individuals in the two cohorts is shown in Table S1. Raw reads in fastq format were fetched from the European Nucleotide Archive (ENA) database using Aspera software. Clinical data of individuals with HCC were obtained from two previous studies by Yoon et al. and Candia et al. For preprocessing, the sequencing adapters were first removed, and the sequencing quality of the reads was checked using fastp. Then the filtered RNA-seq reads were aligned to the human reference genome hg38 using STAR. Raw gene counts were quantified using featureCounts from the Subread toolkit and then transformed into transcripts per kilobase of exon model per million mapped reads (TPM) using an in-house R script.,

GEP and somatic mutation data of the TCGA-LIHC cohort were obtained using the R TCGAbiolinks and TCGAmutations package, respectively., Somatic copy number variation data and clinical information for the TCGA-LIHC cohort were downloaded from the Firebrowse database (http://firebrowse.org/) and University of California, Santa Cruz (UCSC) Xeno (https://xena.ucsc.edu/public).

For drug sensitivity analyses, the raw microarray data (GSE104580) of individuals with HCC who received TACE therapy were analyzed using the R affy package. Additionally, 34 HCC cell lines with gene expression as well as drug sensitivity data (AUC values) were acquired from the Liver Cancer Cell Lines Database (http://lccl.zucmanlab.com).

eRNA quantification

A previously published pipeline, Pipeline for Enhancer Transcription (PET) (http://fun-science.club/PET/), was adopted for eRNA identification and quantification in our study. Briefly, the trimmed fastq reads were aligned to the hg19 reference genome using STAR, followed by quantification of eRNA using featureCounts with an eRNA annotation file from PET. We calculated reads per kilobase of transcript per million mapped reads (RPKM) values from raw count data and filtered out eRNAs whose length was less than 50 bp to avoid mismatches and multi-mapped transcripts, as reported previously. The detected eRNAs were defined as RPKM > 0 in more than 10% of tumor or normal samples and subjected to analyses. The coordinate distribution of eRNAs was illustrated by the R circlize package.

eRNA differential expression analysis and functional annotation

As reported in a previous study, we estimated the difference of eRNA expression at the single and global level. Global eRNA expression was evaluated by aggregating the expression level of single eRNAs in each sample and then averaged by the number of detected eRNAs. For single-eRNA-level analysis, the difference in log2-transformed RPKM value of each eRNA between two groups was evaluated by Wilcoxon test, and the p value was corrected by false discovery rate (FDR) using the Benjamini-Hochberg method for multiple-sample tests. A single eRNA with an absolute log2FC (fold change) value greater than 1 and FDR less than 0.05 was defined as differentially expressed.

To annotate the biological functions of the upregulated eRNAs (log2FC > 1 and FDR < 0.05) in tumors, we first performed Spearman co-expression analysis between eRNAs and PCGs. PCGs with a correlation coefficient greater than 0.6 and p value less than 0.05 were considered positively co-expressed. Then we estimated the enriched KEGG pathways (adjusted p < 0.05) of the co-expressed genes with the R clusterProfiler package.

Identification of eRNA-associated HCC subtypes

NMF subtyping was conducted to subtype the eRNA-related HCC subclasses in both cohorts by the R NMF package. We first screened out eRNAs that were detected in over 50% of individuals and then chose the first 50% of eRNAs (n = 1,244) with a higher median absolute deviation value in the Candia et al. cohort to perform the NMF analysis. The subtyping procedure was also applied to the Yoon et al. cohort with most of the same candidate eRNAs as in the Candia et al. cohort, except for 9 eRNAs that were not detected in the Yoon et al. cohort. The eRNAs used in the NMF analysis are shown in Table S3. In the NMF subtyping procedure with 30 iterations, the optimal number of subtypes was determined by determining the cophenetic score (reproducibility of a model) from 2 to 7 subtypes where the value of the cophenetic score started to fall sharply. PCA was conducted to validate the subtype assignments. Based on the eRNA expression profile, subclass mapping analysis (Submap) was carried out to evaluate the similarity of eRNA subtypes from the two cohorts.

Molecular characterization and immune infiltration estimation of HCC subtypes

Single-sample gene set enrichment analysis (ssGSEA), based on 50 cancer hallmark gene sets and 186 KEGG gene sets, was performed to evaluate the biological properties of HCC subtypes with the R GSVA package. Previously published HCC molecular subtypes were predicted by the MS.liverK algorithm.

To illustrate the immune infiltration of HCC subtypes, we first compared the immune score and stromal score among different HCC subtypes by applying the R estimate package. Then the TME landscape was characterized by inferring the abundance of immune and stromal (endothelial and fibroblast) cells from the transcriptome data via the TIMER2 online tool.74, 75, 76 We calculated the cytolytic activity score using the MCP-counter algorithm.

Prediction of the response to immunotherapy, targeted therapy, and TACE for HCC subtypes

Based on SubMap analysis, we first evaluated the similarity of gene expression patterns between HCC subtypes and previously reported individuals with melanoma who received CTLA-4 blockade or PD-1 blockade immunotherapy to infer the probability of benefits from ICIs. Similarly, the sensitivity to two targeted drugs, sorafenib and cabozantinib, was evaluated according to the data of HCC cell lines. To be more specific, the HCC cell lines were arranged from low to high according to the drug AUC value, and cell lines in the first one-third were regarded as drug sensitive, whereas the last one-third were drug resistant. TACE therapeutic efficacy was also estimated for the identified subtypes based on our microarray data.

Identification of eRNA signature for HCC subtype prediction

For HCC subtype prediction, we first performed differential expression analyses between one subtype and the other two subtypes in every cohort to identify the subtype-specific eRNAs. Only eRNAs with a significant difference (log2FC > 1, FDR < 0.05) were defined as subtype-specific eRNAs. We also identified subtype-specific eRNAs that were shared by two associated subtypes from the two cohorts as the classifier eRNAs. As a result, a 51-eRNA classifier was established. To test the reliability of the classifier, the concordance of predictive subtype, through nearest template prediction (NTP) analysis, was compared with previously identified subtypes based on NMF in the Yoon et al. and Candia et al. cohorts, respectively.

Validation of the eRNA subtypes in the TCGA and SYSUCC HCC cohorts

Using the NTP algorithm, we classified individuals with HCC from the TCGA-LIHC and SYSUCC HCC cohorts into three subtypes based on the 51-eRNA classifier. The eRNA expression profile of the TCGA-LIHC cohort was downloaded from the Cancer eRNA Atlas (https://bioinformatics.mdanderson.org/public-software/tcea). Survival difference, molecular characterization, immune infiltration, and drug sensitivity analyses were also performed as described above.

Correlation of eRNA subtypes with clinical features, CNA burden, gene mutations, and neoantigens

Correlations between clinical features and HCC subtypes were analyzed in the Candia et al. and Yoon et al. cohorts. For gene somatic mutation data analysis in the Candia et al. and TCGA-LIHC cohorts, TMB was calculated with the R Maftools package. The mutation frequency of several HCC-related genes was compared. The neoantigens data for individuals with HCC in the TCGA-LIHC cohort were acquired from a previous study. CNA burden was quantified by CNApp (https://tools.idibaps.org/CNApp/).

IHC

HCC samples were fixed in 4% formaldehyde for 24 h at room temperature, moved into 70% ethanol for 12 h, and then embedded in paraffin. After cutting and baking at 60°C for 60 min for de-paraffinization, the samples were incubated in graded alcohol and then treated for antigen unmasking. For IHC staining, endogenous peroxidases were inactivated by 3% hydrogen peroxide at room temperature (RT) for 15 min. The slides were treated with Tris-EDTA buffer (10 mM, pH 8.0) at 121°C for 10 min for antigen retrieval. Tissues were stained with primary antibodies for 12 h at 4°C. After washing with PBS-Tween, tissues were stained with horseradish peroxidase (HRP)-conjugated secondary antibodies against mouse or rabbit for 30 min at 37°C. Then the samples were co-stained with hematoxylin, followed by washing with PBS. The dilution ratios for IHC were 1:500 for CD8 antibodies (ZSGB-BIO, ZA-0508), 1:500 for ki67 antibodies (Invitrogen, PA5-114437), and 1:500 for HRP-conjugated secondary antibodies (ZSGB-BIO, DS-0003).

Statistical analysis

All statistical analyses and data illustrations were conducted using R software (v.4.0.0). Kaplan-Meier survival analysis and log rank test between different subtypes were performed by utilizing the pairwise_survdiff() function in the R survival package. Cox proportional hazards modeling and ROC curve analyses were conducted using the R survminer and timeROC packages, respectively. Chi-square test or Fisher’s exact test was used to test the difference in categorical data. For differences in continuous data between two or multiple groups, Wilcoxon test or Kruskal-Wallis test was performed. The p values were corrected by FDR using the Benjamini-Hochberg method when multiple comparisons occurred. A two-tailed p value od less than 0.05 was regarded as statistically significant.

Availability of data

The TCGA-LIHC and GEO datasets used in this study are available in public databases. The data of the SYSUCC cohort is available from the corresponding authors upon reasonable request.