Hepatocellular carcinoma (HCC) is one of the most common malignancies worldwide. Costimulatory molecules have been proven to be the foundation of immunotherapy. However, the potential roles of costimulatory molecule genes (CMGs) in HCC remain unclear. Our study is aimed to develop a costimulatory molecule-related gene signature that could evaluate the prognosis of HCC patients.
Based on The Cancer Gene Atlas (TCGA) database, univariate Cox regression analysis was applied in CMGs to identify prognosis-related CMGs. Consensus clustering analysis was performed to stratify HCC patients into different subtypes and compared them in OS. Subsequently, the LASSO Cox regression analysis was performed to construct the CMGs-related prognostic signature and Kaplan–Meier survival curves as well as ROC curve were used to validate the predictive capability. Then we explored the correlations of the risk signature with tumor-infiltrating immune cells, tumor mutation burden (TMB) and response to immunotherapy. The expression levels of prognosis-related CMGs were validated based on qRT-PCR and Human Protein Atlas (HPA) databases.
All HCC patients were classified into two clusters based on 11 CMGs with prognosis values and cluster 2 correlated with a poorer prognosis. Next, a prognostic signature of six CMGs was constructed, which was an independent risk factor for HCC patients. Patients with low-risk score were associated with better prognosis. The correlation analysis showed that the risk signature could predict the infiltration of immune cells and immune status of the immune microenvironment in HCC. The qRT-PCR and immunohistochemical results indicated six CMGs with differential expression in HCC tissues and normal tissues.
In conclusion, our CMGs-related risk signature could be used as a prediction tool in survival assessment and immunotherapy for HCC patients.
Primary liver cancer is an aggressive malignant tumor with high mortality worldwide . Hepatocellular carcinoma (HCC) is the most common histological subtype and the fourth leading cause of cancer-related mortality, accounts for approximately 90% of all primary liver cancer. At present, the traditional treatment methods for HCC are systemic chemotherapy, local ablation, TACE (Transhepatic Arterial Chem Otherapy and Embolization) and surgical resection . However, the therapeutic effect of these methods is away from satisfactory. The effect of anti-tumor agents was not consistent in different clinical trials. Portal vein thrombosis and cholestasis or biliary fistula still needed to be addressed in ablation process, when the target lesion was close to vessels and bile ducts [3, 4]. The recurrence rate of liver cancer after resection was up to 70% at 5 years . In recent years, some clinical trials related immunotherapy showed different outcomes in improving the prognosis of HCC patients [6,7,8]. Therefore, it is urgently required to explore novel prognostic signature for HCC that can predict survival and the response to immunotherapy.
The component of immune microenvironment in HCC is the target for many therapeutic advances, including immunotherapy . Most recently, immunotherapies targeting the adaptive immune system, specifically, T cells, have improved tumor control . Activating T cells involves many signals, among which costimulatory molecules are important [11, 12]. HCC could utilize immune checkpoint and evade anti-tumor immune responses by expressing the corresponding costimulatory ligands . B7-CD28 superfamily is a pivotal signal in co-stimulation of T cell activation, and PD-1/PD-L1 also belong to it, which demonstrated the critical effect of costimulatory molecules in HCC [14, 15]. Besides, accumulating evidence has shown that TNF superfamily, another costimulatory signals, plays a central role in cancer immune regulation . The OX40-OX40L axis, a member of the TNF superfamily, has been shown to improve anti-tumor effects of immune cells and effect for cancer immunotherapy [17,18,19]. Previous studies also have shown that costimulatory molecules can regulate the tumor immune microenvironment (TME), mainly affecting the activation and proliferation of T cells . Thus, these molecules possibly could provide novel insights in TME. However, the functions of costimulatory molecules in HCC remain unclear.
In this systematic study, we evaluated the expression levels of costimulatory molecules genes in HCC tissues and normal tissues from The Cancer Genome Atlas (TCGA) database. Then a costimulatory molecules-related prognostic signature was constructed for HCC patients and we explored the associations between the prognostic signature and clinicopathological features. Furthermore, we also analyzed the potential roles of this prognostic signature in the immune microenvironment, tumor mutation analysis and response to immunotherapy in different subgroups.
Materials and methods
The transcriptomic data and corresponding clinical information of HCC were downloaded from the public The Cancer Genome Atlas (TCGA) data portal (https://portal.gdc.cancer.gov/). A total of 50 normal samples and 374 HCC samples were obtained. Patients with incomplete overall survival (OS) information were excluded. Subsequently, the TCGA cohort was randomly divided into training set (n = 186) and test set (n = 184). There were no significantly differences in clinical characteristics between two sets (Table 1). Furthermore, a total of sixty costimulatory molecules genes (CMGs) were collected from prior reviews [21, 22].
Identification of differentially expressed genes (DEGs)
We utilized “limma” package in R software (version 4.0.4) to identify the differentially expressed genes (DEGs) between all HCC specimens and normal specimens according to the criteria of P-value < 0.05 and |log2 (fold change) |> 1. The DEGs were notated with *** if P < 0.001, ** if P < 0.01 and * if P < 0.05. A PPI network was constructed using the Search Tool for the Retrieval of Interacting Genes (STRING) database (http://www.string-db.org/) to explore the interactions between these DEGs.
Consensus clustering of prognosis-related CMGs
Univariate Cox regression analysis was performed to screen the CMGs with prognostic values in HCC with the cutoff value of P < 0.05. To further elucidate the biological characteristics and prognostic values of CMGs, we employed the “ConsensusClusterPlus” package to cluster the HCC patients into different subgroups . Principal Component Analysis (PCA) was performed using R package to assess the distribution of gene expression among different subtypes. The OS difference between different clusters was verified by the Kaplan–Meier curves. Gene set enrichment analysis (GSEA) was conducted in gene set “h.all.v7.2.symbols.gmt” using Java GSEA software (version 4.1.0) to identify the potential biological processes among different clusters. An enrichment pathway with the normalized P < 0.05 and the false discovery rate (FDR) value < 0.05 were considered as statistically significant.
Construction of costimulatory molecule-related prognostic signature
Patients with HCC were randomly divided into a training set and a test set. The training set was used to construct a prognostic costimulatory molecule-related risk signature of HCC, and the test set and total set were used to validate the prognostic power of this risk signature. The least absolute shrinkage and selection operator (LASSO) penalized Cox proportional hazards regression was performed to narrow down the candidate genes and construct the risk model based on the prognosis-related costimulatory molecule genes using the R package “glmnet” . The penalty parameter (λ) was determined by the minimum criteria. The risk score was calculated with the following formular for each patient: Risk score = expression of gene 1 * coefficient 1 + expression of gene 2 * coefficient 2 + expression of gene 3 * coefficient 3 + … + expression of gene n * coefficient n . Patients were divided into high- and low-risk groups according to the median cutoff of the risk score. The area under the curve (AUC) was calculated between high- and low-risk groups with R package “survivalROC” to validate the prognostic capability. The Kaplan–Meier survival curves of the high- and low-risk groups were plotted using R package “survival” and “survminer” to demonstrate the OS of the patients.
Construction and validation of a nomogram
The nomogram and calibration curves were constructed with R package “rms”. The consistency between the predicted and actual survival of the calibration curves was used to evaluate the accuracy of the nomogram. Meanwhile, the nomogram was examined using the ROC curves.
Functional enrichment analysis
HCC patients were stratified into high- and low-risk groups based on the median risk score. To explore the potential molecular mechanisms of the risk model genes, DEGs between the high- and low-risk groups were identified with the criteria of |log2FC|≥ 1 and FDR < 0.05. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted using the “clusterProfler” package in R software according to the DEGs.
Assessment of immune cell infiltration
The Cell-type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) analysis were performed to estimate the proportions of immune cells infiltration using R package “CIBERSORT” from RNA-sequencing data in TCGA . Wilcoxon rank-sum test was used to examine the differences of infiltrating immune cells in high- and low-risk groups. The tumor microenvironment score was calculated using R package “ESTIMATE” .
The mutation data for HCC patients were downloaded from the TCGA data portal (https://portal.gdc.cancer.gov/). Mutation data were further analyzed using the “maftools” package . We calculated the tumor mutation burden (TMB) score for each HCC patient as follows: (total mutation/total covered bases) × 10^6 .
Immunophenoscore (IPS) could well predict the response of immune checkpoint inhibitors (ICIs). The immunogenicity is determined by four major categories of genes, including effector cells, major histocompatibility complex (MHC) molecules, immunomodulators and immunosuppressive cells. The IPS of a patient can be derived using machine learning without bias. The scores of IPS were calculated using a scale ranging from 0–10 based on representative cell type gene expression z-scores. The IPS of every HCC patient was obtained from The Cancer Immunome Atlas (TCIA) (https://tcia.at/home).
Verification of prognosis-related CMGs expression
Total RNA was extracted from tissue samples using Trizol reagent (Sigma, USA), and then, RNA was reverse transcribed into cDNA with the Evo M-MLV RT Premix (Accurate Biotechnology (Hunan) Co.,Ltd). Quantitative real-time PCR (qRT-PCR) analyses were performed by SYBR Green premix pro Taq HS qRT-PCR kit (Accurate Biotechnology (Hunan) Co.,Ltd) to validate gene expression, and the level of β-Actin served as an internal control. The relative expression was calculated based on the comparative Ct (2−ΔΔCt) method. The primers’ sequences for qRT-PCR are shown in Table 2. The protein expression levels of 6 prognostic gene signatures in normal liver and HCC tissues were determined using the Human Protein Atlas (HPA).
Forty-three matched tumorous and non-tumorous tissue specimens of HCC were collected from The Xijing Hospital of Air Force Medical University during 2017–2018. None of the enrolled patients had received any antitumor agents, such as chemotherapeutic agents, targeted agents, or immunotherapy, prior to surgical resection. The clinicopathological details are shown in Table 3. The research was approved by the Institutional Research Ethics Committees of the Xijing Hospital. Informed consent for publication was obtained from all patients for collection of tissue samples prior to the surgery.
Identification of DEGs between normal and HCC tissues
The flowchart of this study was illustrated in Fig. 1. The expression data of 59 CMGs, including 13 well-defined B7-CD28 family costimulatory molecules and 46 TNF family costimulatory molecules genes, were extracted from The Cancer Genome Atlas (TCGA) database after excluding TNFRSF6B for its low expression. The 59 costimulatory molecule-related genes expression levels were compared between HCC tumor and normal tissues, we identified 40 differentially expressed genes (DEGs) (P < 0.05). Among these DEGs, 11 genes (NGFR, TNFSF11, PDCD1LG2, CD274, TNFRSF1A, TNFRSF11B, TMIGD2, FAS, TNFRSF10D, TNFSF13 and CD86) were down-regulated while 29 genes (TNFRSF17, TNFRSF13B, CD276, TNFRSF12A, LTBR, TNFSF18, EDAR, TNFRSF14, ICOSLG, RELT, CD28, ICOS, LTA, TNFRSF21, TNFRSF10C, VTCN1, TNFRSF11A, LTB, EDA2RC, TLA4, TNFSF9, TNFRSF25, PDCD1, CD70, TNFSF4, TNFRSF9, TNFRSF18, TNFSF15 and TNFRSF4) were up-regulated in tumor tissues (Fig. 2A and Additional file 1: Table S1). The correlation among CMGs were analyzed. The relationships between each two of them were almost positively correlated, TNFRSF13C and TNFRSF13B were most correlated (Cor = 0.93) (Fig. 2B). A protein–protein interaction (PPI) network was performed to further explore the interactions among these CMGs (Fig. 2C). The minimum required interaction score was set at highest confidence 0.9. The result showed that TNF, CD28, CD40, CD80, CTLA4, LTA, TNFRSF10A and TNFSF13B were hub genes.
Consensus clustering of prognosis-related CMGs in HCC
A univariate Cox regression analysis was performed to primary selecting of the survival-related genes from 59 CMGs. A total of 11 CMGs were significantly linked to the prognosis of HCC patients (P < 0.05). Two genes were protective genes with hazard ratio (HR) < 1, while 9 genes were risk factors with HR > 1 among them (Fig. 3A). To explore the associations between the expression of 11 prognosis-related CMGs and HCC subtypes, we performed a consensus clustering analysis to classify all HCC patients based on the expression patterns of 11 prognosis-related CMGs. The empirical cumulative density function (CDF) plot is aimed to determine the optimum cluster number (k) from 2 to 9 for the sample distribution to reach an approximate maximum, which means the maximum stability. When k = 2, the consensus matrix showed that HCC patients could be divided into two non-overlapping subgroups with the highest consensus and the clearest cluster partition (Cluster 1: n = 197, Cluster 2: n = 173) (Fig. 3B and Additional file 3: Figure S1A). The PCA were analyzed to verify the reliability between different subgroups, and Cluster 1 and Cluster 2 could gather together and non-overlapped with each other (Fig. 3C). We compared the OS between two subtypes to better understand the relationships between clustering results and survival outcomes, the Kaplan–Meier curves indicated that Cluster 1 had a better prognosis than Cluster 2 (P = 0.002, Fig. 3D). The clinical features and two clusters were compared with a heatmap. The majority of 11 prognosis-related CMGs had higher expression in Cluster 2. These two clusters were different in grade (P < 0.01), but not with tumor stage, age and gender (Fig. 3E). Furthermore, GSEA analysis showed that oncogenic pathways (apoptosis, G2M-checkpoint, IL2-STAT5-signaling, IL6-JAK-STAT3-signaling, inflammatory response, PI3K-AKT-MTOR-signaling, TNFA-signaling via NF-κB, unfolded protein response) were significantly enriched in Cluster 2 (Additional file 3: Figure S1B).
Construction and verification of costimulatory molecule-related risk signature
To narrow down candidate genes and construct the risk signature, the least absolute shrinkage and selection operator (LASSO) Cox regression analysis was performed in the training set, and 6 of 11 prognosis-related CMGs were identified (Additional file 4: Figure S2). The formula to calculate the risk score as follows: risk score = (0.25584 * TNFSF4) + (−0.29002 * TMIGD2) + (0.13379 * TNFRSF4) + (0.22009 * TNFRSF11B) + (0.40207 * TNFRSF11A) + (−0.78099 * CD40LG). We calculated the risk scores for every HCC patient in the training set according to the above formular. Patients in the training set were divided into high- and low-risk groups based on their median risk sore. A significant difference of OS was observed in different subgroups. High-risk patients had a poorer OS than low-risk groups (P < 0.001) (Fig. 4D). Time-dependent receiver operating characteristic (ROC) analysis was used to evaluate the sensitivity and specificity of the risk signature. The areas under the curve (AUC) were 0.756 at 1-year survival, 0.791 at 3-year survival and 0.729 at 5-year survival (Fig. 4E). We ranked the risk scores of patients and analyzed their distribution in the training set (Fig. 4A). The survival status of HCC patients in the training set was showed on the dot plot (Fig. 4B). The heatmap displayed the expressions of 6 prognosis-related CMGs between two risk groups (Fig. 4C).
To determine the stability of the risk signature, we further verified the predictive capability in the test set and total set. The risk score was calculated for each patient in the test set and total set by the same formular obtained from the training set and the patients were classified into high- and low-risk groups. Similarly, Kaplan–Meier survival curve showed significantly difference in two risk groups among the test set. The OS of the high-risk groups was poorer than that of the low-risk groups (P = 0.019) (Fig. 4I). The 1-year AUC was 0.728, the 3-year AUC was 0.644 and the 5-year AUC was 0.654 (Fig. 4J). The survival status, the distribution of the risk score and the expression heatmap of 6 prognosis-related CMGs in the test set were presented in Fig. 4F–H.
The results in the total set were similar to the training set and test set. Patients in the high-risk group had a significantly shorter prognosis than patients in the low-risk group (P < 0.001) (Fig. 4N). In the total set, the AUC was 0.739 at 1 year, 0.708 at 3 years and 0.662 at 5 years (Fig. 4O). The distribution of the risk score, survival status and the expression patterns of 6 prognosis-related CMGs were showed in Fig. 4K–M.
Independent prognostic value of the risk signature
We performed the univariate and multivariate Cox regression analyses to examine whether the risk score could act as an independent prognosis variable of HCC. Univariate Cox regression analysis showed that pathological tumor stage and risk score were significantly associated with the prognosis (Fig. 5A, C, E). Multivariate Cox regression analyses further identified that the risk score was an independent prognostic factor for OS in the training set, test set and total set (Fig. 5B, D, F).
Correlations between the risk signature and clinicopathological factors
The association between the risk model and clinical characteristics were analyzed to further verify the prognostic value of the risk signature in HCC. The heatmap displayed the expressions of 6 prognosis-related CMGs and the distribution of clinicopathological characteristics in high- and low-risk groups. The risk score was significantly correlated with histological grade, pathological T stage and clinical stage (Fig. 6A). Differences in clinicopathological factors between high- and low-risk groups were showed in Fig. 6B. The risk score of patients with stage III-IV was higher than that of patients with stage I-II (P = 0.0014). Patients with G3–4 showed a remarkably higher risk score than those with G1–2 (P < 0.001). The risk score in T3–4 patients was significantly higher than that observed in T1–2 patients (P = 0.0015). Nonetheless, there were no significant differences among age, gender, N stage and M stage. Furthermore, patients were divided into different subgroups according to the following clinical variables, including age (≤ 65 and > 65), gender (female and male), clinical stage (stage I-II and stage III-IV), grade (G1–2 and G3–4), T stage (T1–2 and T3–4), N stage (N0 and N1–3) and M stage (M0 and M1). The correlation between the risk score and clinicopathological features on OS was explored. Survival analysis manifested that higher risk score were correlated with poor prognosis in age (age ≤ 65 with P = 0.005 and age > 65 with P < 0.001), gender (P < 0.001 in male), stage (stage I-II with P = 0.002 and stage III-IV with P = 0.026), grade (P < 0.001 in G1–2 and P = 0.029 in G3–4), T stage (P < 0.001 in T1–2 and P = 0.01 in T3–4), stage N0 (P < 0.001), and stage M0 (P < 0.001) (Fig. 6C).
Construction of a novel nomogram
We constructed a nomogram to predict the survival rates for HCC patients based on age, gender, histological grade, TNM stage, tumor stage and risk score (Fig. 7A). The calibration curves of the nomogram indicated good consistency between the predicted survival rate and actual 1-, 3- and 5-year survival rate (Fig. 7B). The AUCs of risk score and tumor stage were 0.739 and 0.671 in 1-year, 0.698 and 0.680 in 3-year, 0.638 and 0.663 in 5-year, respectively (Fig. 7C). These findings suggested that the risk signature might be reliable to predict the OS for HCC patients.
Functional enrichment analyses based on the risk signature
To explore the potential biological processes for the prognostic risk signature, a total of 474 DEGs were obtained in high- and low-risk groups with the criteria FDR < 0.05 and |log2FC |≥ 1. Among them, 63 genes were downregulated in high-risk group, while 411 genes were upregulated (Additional file 2: Table S2). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analyses were carried out based on the DEGs. The results of GO analysis indicated that the DEGs were mainly related to nuclear division. KEGG analysis showed that DEGs were mostly enriched in cell cycle (Additional file 5: Figure S3).
Association between risk signature and tumor immune microenvironment
The differences of tumor-infiltrating immune cells between high- and low-risk groups were analyzed to explore the correlations between the prognostic risk signature and tumor immune microenvironment (TIME). Additional file 6: Figure S4 displayed the abundance of 22 immune cells between high- and low-risk subgroups. Among 22 immune cell types, memory B cells and macrophage M0 were positively correlated with the risk score, while the abundance of naïve B cells, plasma cells and regulatory T cells (Tregs) were significantly enriched in low-risk group (Fig. 8A). Furthermore, the relative proportion of naïve B cells, resting memory CD4 T cells, activated memory CD4 T cells, regulatory T cells, gamma delta T cells, macrophage M1 and resting mast cells were significantly associated with OS (Fig. 8B).
To further explore the relationship between the risk signature and immune status, we performed the expression profiles of 29 immune signature gens sets (16 types of immune cells and 13 immune-related pathways) in high- and low-risk groups using the single-sample gene set enrichment analysis (ssGSEA). The heatmap dis1layed the significant differences in immune status between high- and low-risk samples (Fig. 9A). The low-risk subgroup showed higher levels of infiltration of immune cells and higher activity of immune-related pathways (Fig. 9B). We found that the immune score and stromal score were higher in low-risk groups, while the tumor purity was significantly lower in low-risk subgroup (Fig. 9C).
Differences in molecular characteristics between high- and low-risk groups
We evaluated the relationship between mutation characteristics and the risk signature in TCGA HCC patients with available somatic mutation data. TMB was higher in high-risk patients in spite of no significant difference (Fig. 10B). We also identified the top 20 genes with the highest mutation rates in high- and low-risk subgroups (Fig. 10A). Additionally, we explored the association between immunophenoscore (IPS) and risk signature to predict the potential clinical efficacy and the response to ICI therapy in HCC patient. The IPS, IPS-CTLA4, IPS-PD1-PD-L1- PD-L2, and IPS-PD1-PD-L1-PD-L2-CTLA4 blocker were significantly higher in low-risk group, implying that HCC patients with low-risk score could benefit more from ICI therapy than high-risk patients (Fig. 10C).
Verification of prognostic CMGs
We verified the expression of the CMGs (TNFSF4, TNFRSF4, TMIGD2, TNFRSF11A, TNFRSF11B, CD40LG) in 43 pairs of tumorous and non-tumorous tissue specimens from patients with HCC using qRT-PCR analysis (Fig. 11). The results of qRT-PCR showed that the expression of TNFSF4, TNFRSF4, TNFRSF11A and CD40LG was higher in HCC tissues compared to normal tissues. However, the mRNA expression of TNFRSF11B and TMIGD2 was higher in normal tissues, which was consistent with the results of bioinformatic analysis. The protein expression of the 6 prognostic gene signatures in HCC tissues and normal liver tissues was verified in HPA online database (Fig. 12). The results showed that the expression of TNFRSF11A was increased and the expression of TNFRSF11B was decreased in HCC tissues. TNFRSF4 and TMIGD2 were not detected in hepatocytes. In addition, the staining results of TNFSF4 and CD40LG did not reach significant difference according to HPA database.
Preliminary data from trials of ICIs in the treatment of HCC led to encouraging results. Nevertheless, with the rapid augment in the utilization of ICIs, immune-related adverse events of HCC arose [30, 31]. Multiple studies found that the usage of ICIs strategies targeting costimulatory molecules for HCC management was promising . Therefore, it was necessary to improve the effect on immunotherapy by selecting the suitable HCC patients according to costimulatory molecules expression patterns. In this study, we analyzed the mRNA expression patterns of costimulatory molecules-related in HCC and selected six genes with prognostic values. Then, we constructed the first costimulatory molecule-related prognostic signature for HCC patients. We found that prognostic signature was strongly associated with clinical characteristics. Additionally, our signature was significantly correlated with tumor immune microenvironment and the response to immunotherapy. Univariate and multivariate cox regression analysis indicated our signature could be an independent prognostic factors for the survival of HCC patients. These findings suggested that CMGs risk signature may indicate some insights to personalized targeted treatment in clinical practice.
Costimulatory molecules played an important role in immunotherapy . Recent findings demonstrated that CD28 co-stimulation was necessary for responses to PD-1 blockade in tumor rejection . Thus, understanding the states of costimulatory molecules in HCC patients will help us determine which patients might benefit in immunotherapy. To explore the expression levels of costimulatory molecules in HCC, we acquired 13 members of the B7-CD28 family and 46 members of the TNF family for HCC patients. Six costimulatory molecular genes (TNFSF4, TNFRSF4, TMIGD2, TNFRSF11A, TNFRSF11B, CD40LG) with prognostic values were selected. The TNFRSF4-TNFSF4 pathway provided crucial co-stimulatory signals for CD4+T cell responses . Previous study showed that TNFSF4 was closely related to the unfavorable prognosis of HCC patients . In addition, TNFRSF4 was overexpressed in HCC, associated with a more aggressive phenotype and the activation of multiple immunosuppressive pathways . A phase I clinical research also supported that Ivuxolimab (a TNFRSF4 agonist) showed well tolerance and effective anti-tumor capacity in locally advanced or metastatic cancers, including HCC . Consequently, treatment targeting TNFRSF4-TNFSF4 should be considered in the future. TMIGD2 was mainly expressed in tissue-resident lymphocyte T cells, related to improved tumor prognosis . The different interaction between TMIGD2 and B7-H5 have been identified in certain cancers, such as lung cancer, osteosarcoma, oral squamous cell carcinoma (OSCC), colorectal cancer (CRC) and glioma . The TMIGD2-B7-H5 interaction was involving in Akt-dependent signal pathway, which was a recognized regulation in tumor . Previous study reckoned that HCC patients were likely to be benefited from Everolimus (involved in the Akt-MTOR pathway) only after molecular screening . As such, we could make a reasonable speculation that molecular selection was very necessary for the individualized treatment of HCC and our prognostic signature might provide a new method in distinguishing suitable patients. Of note, our results firstly revealed that TMIGD2 was highly expressed in HCC with favorable prognosis. TNFRSF11A, also known as RANK, was significantly up-regulated in HCC, and can lead directly to migration and invasion by its ligand . Interestingly, genetic deletion of TNFRSF11A in thymic epithelial cells resulted in impaired thymic involution and blunted expansion of natural regulatory T (Treg) cells . Additionally, study showed that HCC patients with high serum TNFRSF11B, also known as osteoprotegerin (OPG), level had poorer survival rates compared with HCC patients with low OPG level . CD40 ligand-expressing dendritic cells could induce regression of HCC . In addition, some clinical trials targeting the agonist or antagonist of CD40/CD40LG has showed promising results in different malignancies. Study showed that adenoviral vectors expressing CD40 ligand (AdCD40L) were safe in vivo and could reduce the number of tumor cell infiltration in bladder cancer . Besides, AdCD40L intratumoral injection increased T-effector/T-regulatory cell ratio by improving systemic immune condition, which was related favorable survival in malignant melanoma patients . Thus, combined other scholar studies and our own bioanalysis, it was possible to improve the prognosis of HCC patients with a similar approach. Moreover, the expression levels of six prognostic genes were verified using qRT-PCR and immunohistochemistry. However, the protein expression of prognostic genes signatures was not completely similar with our previous results, may partly owing to different race and clinical characteristics. At the same time, it also explained that why targeted therapy and immunotherapy did not work for all HCC patients and tumor heterogeneity should be considered in treatment practice. With these six costimulatory molecular genes elucidated in immunity, we hope that the signature constructed by these could predict the response to immunotherapy for HCC patients.
Tumors were complex ecosystems, defined by spatiotemporal interactions among heterogeneous cell types . Subsequently, we compared the associations between our signature and tumor immune microenvironment, the immune cell infiltration and tumor mutation profiles in high-risk and low-risk patients. Our results showed that naïve B cells, plasma cells and regulatory T cells (Tregs) were significantly enriched in low-risk groups. Much of researches regarded Tregs as an immunosuppressive cell, posing anti-tumor immunity in various cancers . Nevertheless, some scholars stated that inhibiting the expression of PD-1 promoted other immune checkpoints, resulting in impaired immune killing ability [50, 51]. Accordingly, fewer Tregs cells in HCC patients with poor prognosis, indicated that those were more likely to be inhibited by PD-1 and activated more immune checkpoints. Correspondingly, high-risk subgroup manifested lower levels of infiltration of immune cells, implicating less process in immune activation.
Tumor mutation burden (TMB) is emerging as a potential biomarker, and participated in immunotherapy-related pathway [52, 53]. We found that the tumor mutation burden (TMB) in the high-risk group was higher than that in the low-risk group with no significant, partly due to the small sample size. TP53 mutation frequency was evidently higher in high-risk group (frequency rate 42%) than low-risk group (frequency rate 14%) according to our mutation results, suggesting more increases genomic instability and complicated major pathway signaling changes in HCC. Additionally, it was important to highlight that different microenvironment-based immune subtypes, based on gene profiling or signatures, and other molecular features, may help identify subgroups of patients more likely to benefit from specific therapies . Some scholars have found that immune-excluded tumors in HCC were proposed to be primarily resistant to ICIs . IPS could predict the response to immunotherapy in cervical cancer and HCC. The prediction of IPS has been demonstrated in different studies [56, 57]. In the present study, we found low-risk group tended to have higher IPS-CTLA4, IPS-PD1/PD-L1/ PD-L2, and IPS-PD1/PD-L1/PD-L2+CTLA4, implying that HCC patients with low-risk score could benefit more from immunotherapy than high-risk patients. Therefore, our signature was of great help to clinical immunotherapy decision.
However, there were some limitations in this study. Firstly, we did not explore the exact function of six costimulatory molecule genes in HCC. Thus, it was still necessary to clarify the mechanism of them in the future. Secondly, it was inevitable that there were limited clinical information for HCC patients in public datasets, so the values of the prognostic signature needed to be determined by experimental and prospective studies. Moreover, the risk signature for evaluating the response to immunotherapy was restricted to costimulatory molecule genes and tumor immune microenvironment was highly heterogeneous. Therefore, the prognostic information for HCC patients with immunotherapy were needed to validate the prediction power of our signature clinically.
In our study, we first elucidated the expression of costimulatory molecules for HCC patients, and constructed a six CMGs prognostic signature. The costimulatory molecular-related signature could stratify patients into different subsets with adverse clinical outcomes. In addition, immunotherapy response prediction by our signature explained disparate effect on HCC patients. Consequently, we believed our research manifested the capacity of costimulatory molecules and provided clinicians with applicable treatment.
Availability of data and materials
Gene expression profiles, clinical information and mutation data of LIHC in this study are available from the public database (TCGA, https:// portal. gdc. cancer. gov/).
Sung H, Ferlay J, Siegel RL, Global cancer statistics, et al. GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2020. https://doi.org/10.3322/caac.21660.
Ding J, Jing X, Liu J, et al. Complications of thermal ablation of hepatic tumours: comparison of radiofrequency and microwave ablative techniques. Clin Radiol. 2013. https://doi.org/10.1016/j.crad.2012.12.008.
Franssen B, Alshebeeb K, Tabrizian P, et al. Differences in surgical outcomes between hepatitis B- and hepatitis C-related hepatocellular carcinoma: a retrospective analysis of a single North American center. Ann Surg. 2014. https://doi.org/10.1097/SLA.0000000000000917.
Zhu AX, Finn RS, Edeline J, et al. Pembrolizumab in patients with advanced hepatocellular carcinoma previously treated with sorafenib (KEYNOTE-224): a non-randomised, open-label phase 2 trial. Lancet Oncol. 2018. https://doi.org/10.1016/s1470-2045(18)30351-6.
Lee MS, Ryoo B-Y, Hsu C-H, et al. Atezolizumab with or without bevacizumab in unresectable hepatocellular carcinoma (GO30140): an open-label, multicentre, phase 1b study. Lancet Oncol. 2020. https://doi.org/10.1016/s1470-2045(20)30156-x.
Lu Y, Zhang M, Wang S, et al. p38 MAPK-inhibited dendritic cells induce superior antitumour immune responses and overcome regulatory T-cell-mediated immunosuppression. Nat Commun. 2014. https://doi.org/10.1038/ncomms5229.
Aberg KM, Radek KA, Choi EH, et al. Psychological stress downregulates epidermal antimicrobial peptide expression and increases severity of cutaneous infections in mice. J Clin Invest. 2007. https://doi.org/10.1172/JCI31726.
Shibahara I, Saito R, Zhang R, et al. OX40 ligand expressed in glioblastoma modulates adaptive immunity depending on the microenvironment: a clue for successful immunotherapy. Mol Cancer. 2015. https://doi.org/10.1186/s12943-015-0307-3.
Aye L, Song X, Yang J, et al. Identification of a costimulatory molecule gene signature to predict survival and immunotherapy response in head and neck squamous cell carcinoma. Front Cell Dev Biol. 2021. https://doi.org/10.3389/fcell.2021.695533.
Hua X, Ge S, Zhang J, et al. A costimulatory molecule-related signature in regard to evaluation of prognosis and immune features for clear cell renal cell carcinoma. Cell Death Discov. 2021. https://doi.org/10.1038/s41420-021-00646-2.
Becht E, Giraldo NA, Lacroix L, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 2016. https://doi.org/10.1186/s13059-016-1070-5.
Finn RS, Ryoo BY, Merle P, Bouattour M, Lim HY, Breder V, Edeline J, Chao Y, Ogasawara S, Yau T, Garrido M. Pembrolizumab as second-line therapy in patients with advanced hepatocellular carcinoma in KEYNOTE-240: a randomized, double-blind, phase III trial. J Clin Oncol. 2020. https://doi.org/10.1200/JCO.19.01307.
Chinnadurai R, Scandolara R, Alese OB, et al. Correlation patterns among B7 family ligands and tryptophan degrading enzymes in hepatocellular carcinoma. Front Oncol. 2020. https://doi.org/10.3389/fonc.2020.01632.
Xie K, Xu L, Wu H, et al. OX40 expression in hepatocellular carcinoma is associated with a distinct immune microenvironment, specific mutation signature, and poor prognosis. Oncoimmunology. 2018. https://doi.org/10.1080/2162402X.2017.1404214.
Diab A, Hamid O, Thompson JA, et al. A Phase I, open-label, dose-escalation study of the OX40 agonist ivuxolimab in patients with locally advanced or metastatic cancers. Clin Cancer Res. 2021. https://doi.org/10.1158/1078-0432.CCR-21-0845.
Zhong C, Lang Q, Yu J, et al. Phenotypical and potential functional characteristics of different immune cells expressing CD28H/B7-H5 and their relationship with cancer prognosis. Clin Exp Immunol. 2020. https://doi.org/10.1111/cei.13413.
Koeberle D, Dufour JF, Demeter G, et al. Sorafenib with or without everolimus in patients with advanced hepatocellular carcinoma (HCC): a randomized multicenter, multinational phase II trial (SAKK 77/08 and SASL 29). Ann Oncol. 2016. https://doi.org/10.1093/annonc/mdw054.
Song FN, Duan M, Liu LZ, et al. RANKL promotes migration and invasion of hepatocellular carcinoma cells via NF-kappaB-mediated epithelial-mesenchymal transition. PLoS ONE. 2014. https://doi.org/10.1371/journal.pone.0108507.
Zhang C, Lin J, Ni X, et al. Prognostic value of serum osteoprotegerin level in patients with hepatocellular carcinoma following surgical resection. Front Oncol. 2021. https://doi.org/10.3389/fonc.2021.731989.
Gonzalez-Carmona MA, Lukacs-Kornek V, Timmerman A, et al. CD40ligand-expressing dendritic cells induce regression of hepatocellular carcinoma by activating innate and acquired immunity in vivo. Hepatology. 2008. https://doi.org/10.1002/hep.22296.
Schiza A, Wenthe J, Mangsbo S, et al. Adenovirus-mediated CD40L gene transfer increases Teffector/Tregulatory cell ratio and upregulates death receptors in metastatic melanoma patients. J Transl Med. 2017. https://doi.org/10.1186/s12967-017-1182-z.
Ellestad KK, Thangavelu G, Ewen CL, et al. PD-1 is not required for natural or peripherally induced regulatory T cells: Severe autoimmunity despite normal production of regulatory T cells. Eur J Immunol. 2014. https://doi.org/10.1002/eji.201444688.
Huang RY, Francois A, McGray AR, et al. Compensatory upregulation of PD-1, LAG-3, and CTLA-4 limits the efficacy of single-agent checkpoint blockade in metastatic ovarian cancer. Oncoimmunology. 2017. https://doi.org/10.1080/2162402X.2016.1249561.
Huang T, Yan T, Chen G, et al. Development and validation of a gene mutation-associated nomogram for hepatocellular carcinoma patients from four countries. Front Genet. 2021. https://doi.org/10.3389/fgene.2021.714639.
Sia D, Jiao Y, Martinez-Quetglas I, Kuchuk O, Villacorta-Martin C, de Moura MC, Putra J, Camprecios G, Bassaganyas L, Akers N, Losic B. Identification of an immune-specific class of hepatocellular carcinoma, based on molecular features. Gastroenterology. 2017. https://doi.org/10.1053/j.gastro.2017.06.007.
Galarreta M, Bresnahan E, Molina-Sanchez P, et al. Beta-catenin activation promotes immune escape and resistance to anti-PD-1 therapy in hepatocellular carcinoma. Cancer Discov. 2019. https://doi.org/10.1158/2159-8290.CD-19-0074.
Yang S, Wu Y, Deng Y, et al. Identification of a prognostic immune signature for cervical cancer to predict survival and response to immune checkpoint inhibitors. Oncoimmunology. 2019. https://doi.org/10.1080/2162402X.2019.1659094.
This study was funded by the National Key Research and Development Project (Nos. 2020YFA0710803), Natural Science Foundation of China (Nos. 81770569, 81870421, and 81773072) and the International Cooperation and Exchange of the National Natural Science Foundation of China (No. 81820108005).
Yinan Hu and Jingyi Liu contributed equally to this article
Authors and Affiliations
Institute of Digestive Diseases, Xijing Hospital, Air Force Medical University, Xi’an, 710032, Shaanxi, China
Yinan Hu, Jiahao Yu, Fangfang Yang, Miao Zhang, Yansheng Liu, Shuoyi Ma, Xia Zhou, Jingbo Wang & Ying Han
Department of Radiation Oncology, Xijing Hospital, Air Force Medical University, Xi’an, 710032, Shaanxi, China
Study concept and design (YH), acquisition of data (YNH, JYL), analysis and interpretation of data (YNH, JHY, JYL), drafting of the manuscript (YNH, JHY), critical revision of the manuscript for important intellectual content (MZ, YSL, JBW, YH), administrative, technical, or material support, study supervision (XZ, FFY, SYM). All authors read and approved the final manuscript.
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 conformed to the provisions of the Declaration of Helsinki (as revised in 2013). The research was approved by the Institutional Research Ethics Committees of the Xijing Hospital. Informed consent for publication was obtained from all patients for collection of tissue samples prior to the surgery.
Consent for publication
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
: Consensus clustering based on the CMGs. (A) Cumulative distribution function (CDF) curves in consensus clustering and relative changes in area under CDF curves for k = 2 to 10. (B) The Gene Set Enrichment Analysis of the oncogenic pathways in Cluster 2.
: The abundance of 22 immune cells between the high- and low-risk groups.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
Hu, Y., Liu, J., Yu, J. et al. Identification and validation a costimulatory molecule gene signature to predict the prognosis and immunotherapy response for hepatocellular carcinoma.
Cancer Cell Int22, 97 (2022). https://doi.org/10.1186/s12935-022-02514-0