Yuhong Tang,Jun Cao,Rui Peng,Xingkang Mao,Bingbing Su,Hao Tang,Daoyuan Tu,Jie Zhou,Guoqing Jiang,Shengjie Jin,Qian Wang,Chen Zhang,Renjie Liu,Chi Zhang,Dousheng Bai
Screening and Verification of Key Ubiquitination Genes Related to Immune Infiltration in Stage III/IV Hepatocellular Carcinoma Tang et al
ubiquitination,growth factor receptor-binding protein 2,liver cancer,tumour-infiltrating immune cells
Introduction Immune checkpoint therapy (ICIs) effectively improves the prognosis of advanced (stage III/IV) hepatocellular carcinoma (HCC) patients. However, its objective response rate (ORR) is below 20%, significantly limiting ICI use in advanced HCC patients. The level of tumour immune infiltration influences ICI response rate. Recent studies have found ubiquitinase to be an important factor that regulates tumour immune infiltration. Therefore, the aim of this study is to explore the key ubiquitination genes that regulate immune infiltration in advanced HCC and further validate them. Methods A biotechnological process was performed as a means of classifying 90 advanced HCC patients into three immune subtypes and identifying associations with immune infiltration in the co-expressed modules. Ubiquitination-related genes were then screened with WGCNA. Gene enrichment analysis was performed for the target module and 30 hub genes were screened out by protein–protein interaction network (PPI). ssGSEA, single-gene sequencing and the MCP counter were used for exploring immune infiltration. TIDE score was applied for predicting drug efficacy and GSEA was used for exploring potential pathways. Finally, GRB2 expression in HCC tissue was validated by in vitro experiments. Results GRB2 expression was found to have a significant correlation with the pathological stage and prognosis of HCC patients and a positive correlation with immune infiltration and tumour mutation burden (TMB). In addition, significant correlations with the efficacy of ICIs, sorafenib and transarterial chemoembolization (TACE) were identified. GRB2 was found to be most significantly associated with the JAK-STAT signalling pathway and cytosolic DNA sensing pathway. Finally, it was found that GRB2 expression is closely related to the prognosis, tumour size and TMN stage. Conclusion A significant association was observed between the ubiquitinated gene GRB2 and the prognosis and immune infiltration of advanced HCC patients and it may potentially be used for predicting therapy efficacy in advanced HCC patients in the future.
Screening and Verification of Key Ubiquitination Genes Related to Immune Infiltration in Stage III/IV Hepatocellular Carcinoma Tang et al

Immune checkpoint therapy (ICIs) effectively improves the prognosis of advanced (stage III/IV) hepatocellular carcinoma (HCC) patients. However, its objective response rate (ORR) is below 20%, significantly limiting ICI use in advanced HCC patients. The level of tumour immune infiltration influences ICI response rate. Recent studies have found ubiquitinase to be an important factor that regulates tumour immune infiltration. Therefore, the aim of this study is to explore the key ubiquitination genes that regulate immune infiltration in advanced HCC and further validate them.

A biotechnological process was performed as a means of classifying 90 advanced HCC patients into three immune subtypes and identifying associations with immune infiltration in the co-expressed modules. Ubiquitination-related genes were then screened with WGCNA. Gene enrichment analysis was performed for the target module and 30 hub genes were screened out by protein–protein interaction network (PPI). ssGSEA, single-gene sequencing and the MCP counter were used for exploring immune infiltration. TIDE score was applied for predicting drug efficacy and GSEA was used for exploring potential pathways. Finally, GRB2 expression in HCC tissue was validated by in vitro experiments.

GRB2 expression was found to have a significant correlation with the pathological stage and prognosis of HCC patients and a positive correlation with immune infiltration and tumour mutation burden (TMB). In addition, significant correlations with the efficacy of ICIs, sorafenib and transarterial chemoembolization (TACE) were identified. GRB2 was found to be most significantly associated with the JAK-STAT signalling pathway and cytosolic DNA sensing pathway. Finally, it was found that GRB2 expression is closely related to the prognosis, tumour size and TMN stage.

A significant association was observed between the ubiquitinated gene GRB2 and the prognosis and immune infiltration of advanced HCC patients and it may potentially be used for predicting therapy efficacy in advanced HCC patients in the future.

Liver cancer is the third most common cause of cancer death globally with hepatocellular carcinoma (HCC) accounting for the majority of liver cancer.1,2 Early HCC (stage I/II) treatment mainly includes radical surgical resection, radiofrequency ablation and liver transplantation.3,4 However, there is a lack of typical clinical symptoms with early-stage HCC and over 60% of HCC patients are initially diagnosed at an advanced stage, meaning that the best opportunity for surgery is lost and resulting in a low five-year survival rate (below 16%).5 According to the American Gastroenterological Association (AGA) Clinical Practice Guideline on Systemic Therapy for HCC, patients in stage III or IV were only previously able to receive conservative treatment, including molecular targeted drugs (sorafenib, etc.), TACE and radiation therapy.6 However, these therapies demonstrated limited effect in improving advanced HCC patient prognosis. In recent years, immune checkpoint inhibitors (ICIs) that target PD-1, PD-L1 and CTLA-4 pathways have proven to be effective in many advanced cancers.7,8 Recent studies have reported that a combination of camrelizumab and apatinib has outstanding efficacy in the treatment of advanced liver cancer, with tumour ORR of 34.3% (95% CI: 23.3–46.6%) and a one-year survival rate of 74.7% (95% CI, 62.5–83.5%).9,10 Furthermore, a combination of tremelimumab and TACE treatment was found to increase T lymphocyte accumulation in distant untreated lesions.11 Treatment combining oral hydrogel and ICIs can decrease tumour growth by reversing the immunosuppressive tumour microenvironment (TME).12 Therefore, multiple immune cell subsets infiltrate the TME and play an important role in the response rate of ICIs. It has also been found that ubiquitinated proteases play a key regulatory role in tumour immune infiltration, affecting ICI efficacy.12 Therefore, screening the key molecules of ubiquitinated proteases that regulate immune infiltration in advanced HCC may provide new targets for improving ICI efficacy among HCC patients. Ubiquitinated proteases (E1, E2, E3) are enzymes that catalyse protein degradation by linking ubiquitin as a means of targeting proteins in a three-step cascade and are closely related to antitumour immunity.13,14 Increasing evidence that ubiquitinated proteases play a key role in TME has been identified in recent years. Yu et al reported that ubiquitinase RNF126 affects the TME and promotes nasopharyngeal carcinoma progression through the regulation of PTEN ubiquitination.15 Wang et al discovered that an overexpression of USP24 in M2 macrophages can promote IL-6 transcription and tumour malignancy.16 It was previously predicted that the ubiquitination function gene CDC20 is related to HCC prognosis and its expression has a positive correlation with the infiltration of B cells, CD8 + T cells, CD4 + T cells and DCs by bioinformatics methods.17 The important role of RNF128 in HCC was also identified and it was found that RNF38 promotes TGF-β signalling in HCC by ubiquitinating and degrading AHNAK.18,19 Recently, immune infiltration was used for constructing a patient prognostic model of the ubiquitinates TRIM family as a means of predicting effective prognostic markers and therapeutic targets for HCC.20 These studies suggest that ubiquitinated proteases play a crucial role in tumourigenesis and immune infiltration. Therefore, an exploration of the key molecules of ubiquitinated genes that regulate the immune infiltration of HCC can provide new targets for improving ICI efficacy among HCC patients. In this study, the immune infiltration-related ubiquitinase GRB2 was screened using various bioinformatics approaches. Subsequently, the high expression of GRB2 was found to be associated with a worse prognosis and higher immune cell infiltration levels. In addition, GRB2 expression was found to be closely related to sorafenib, TACE and ICI treatment efficiency among HCC patients. Generally, GRB2 could be a novel predictive biomarker for the prediction of therapy efficacy among advanced HCC patients in the future.

A total of 1597 ubiquitin-protein ligases were obtained from the Integrated Annotations for Ubiquitin and Ubiquitin-like Conjugation Database (IUUCD) (http://iuucd.biocuckoo.org/species.php?spe=Homo_sapiens).21 Gene expression data and the corresponding clinical information of 90 advanced HCC patients were obtained from the International Cancer Genome Consortium (ICGC) database (https://dcc.icgc.org/projects/LIRI-JP) and TCGA database (https://tcga-data.nci.nih.gov/tcga/). In addition, single-cell RNA sequencing data was taken from the “LIHC GSE140228 10X” and “LIHC GSE98638” datasets of the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo).

One hundred and forty pairs of HCC tissues were collected from Yangzhou Clinical Medical College. Thirty pairs were used for quantitative real-time polymerase chain reaction (qRT-PCR) and eight pairs were chosen for Western blotting. All patients in the immunohistochemical staining chip received a radical resection of liver cancer in the Department of Hepatobiliary Surgery, Clinical Medical College, Yangzhou University and the clinical data was complete. Two clinical pathologists independently confirmed the HCC diagnosis. Clinicopathological information was obtained for correlation analysis. All HCC subjects gave their written informed consent in accordance with the Declaration of Helsinki and studies were approved by the Institutional Review Board at Clinical Medical College, Yangzhou University.

Single-sample gene set enrichment analyses (ssGSEA) derived from the gene set variation analysis (GSVA) package was used for scoring immune infiltration in 90 advanced HCC. In addition, the cumulative distribution function (CDF) method was used to determine the optimal grouping according to the CDF delta region. In order to verify differences in the immunity of clusters, the ESTIMATE package for R (https://sourceforge.net/projects/estimateproject) was used for calculating the tumour purity, stromal score, immune score and ESTIMATE score of each HCC tumour sample. The CIBERSORT package was then used for assessing the distribution of 24 immune cell types in each sample. Results with p < 0.05 were then used in further analysis. The R package WGCNA was used for co-expression analysis. Based on the above grouping, selected genes relating to ubiquitination were constructed in a sample clustering tree. Differential genes were then screened as a means of constructing a matrix of similarity by Pearson correlation analysis. A weighted adjacency matrix was then constructed based on soft threshold β. Finally, the 1597 ubiquitination-related genes were divided into four modules and the module with the highest correlation with immune infiltration was identified. The clusterProfiler package was used for performing enrichment analysis of the selected module, based on Gene Ontology (GO) function annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG). On the STRING, module genes in the selected module were utilised for the construction of a PPI network. An interaction score >0.9 was taken as the threshold and the degree score for each gene was determined by Network Analyzer in Cytoscape (https://cytoscape.org).

The correlation between GRB2 expression and immune infiltration cells was verified using ssGSEA and MCP counter was used to score immune infiltration in the TCGA and ICGC databases. The R package heatmap was then used for visualising their correlations according to the coefficient. The data of mutations from the TCGA database were downloaded and visualised using the R maftools package and a p-value of 0.05 was set as the statistical significance threshold.

The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm22 was used to predict the response efficiency to anti-PD-1. “GSE109211” and “GSE104580” from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo) were used for evaluating sorafenib and TACE efficacy. All HCC samples in the TCGA and ICGC databases were divided into a GRB2 low group and GRB2 high group based on the median expressions GRB2. Single-gene GSEA was used for exploring the significant biological pathways between the two groups through GSEA_4.0.1 software (https://www.gsea-msigdb.org/gsea/index.jsp). All results were plotted using the ggplot2 package.

qRT-PCR assay, immunohistochemistry (IHC) and Western blot methods were conducted based on previous literature.22 The following primers were used: GAPDH forward 5′-GT CTCCTCTGACTTCAACAGCG-3′ and reverse 5′-ACCACCCTGTTGCTGTAGCCAA3′; GRB2 forward 5′-TCCTCTGGGTGGTGAAGTTCAATTC-3′ and reverse 5′-GCTGTGGCACCTGTTCTATGTCC3′, and the primary antibodies anti-Tubulin (1:5000) and anti-GRB2 (1:1000) were obtained from Proteintech (Wuhan, China).

All statistical analytical methods were based on IBM SPSS statistical software (version 22.0) and R software (4.0.5). Spearman correlation was used for determining the association between GRB2 expression and immune infiltration cells. In this study, the scatter plot of qRT-PCR results was performed using GraphPad Prism 8 software. A p-value < 0.05 was considered to be statistically significant.

Firstly, ssGSEA was used for scoring immune infiltration in 90 advanced HCC patients from the TCGA database. The optimal number of clusters was determined by the consistent cumulative distribution function graph and the delta area plot. The final partitioning result with K = 3 clusters was obtained (Figure 1A and B). The hierarchical clustering showed three clusters with significant interconnectivity (Figure 1C). In order to determine the relationship of immune scores between clusters in the TCGA database, the ESTIMATE package was used for evaluating immune scores for three clusters and cluster 1 was found to be the immune desert type, cluster 2 the intermediate type and cluster 3 the immune enrichment type (Figure 2A). To validate the above results, it was found that most human leukocyte antigen (HLA) genes exhibited significantly higher expression levels in Immunity_H but lower expression levels in Immunity_L (p < 0.05) (Figure 2B). In addition, the CIBERSORT algorithm quantified the immune cell infiltration level. The results showed the levels of plasma cells, CD8+T cells, memory resting CD4+T cells, and active NK cells to be remarkably increased in cluster 3 (p < 0.001) (Figure 2C). The results suggest that the grouping was representative and advanced HCC patients were further analysed based on these clusters.

A total of 1597 ubiquitination-related genes were first obtained from the IUUCD database and a WGCNA network was constructed. β = 5 was then chosen as the optimal soft-thresholding parameter (Figure 3A). A module–immune cell correlation heat map was drawn to explore the correlation between the four modules and immune cell infiltration. Four modules with co-expressed genes were found to be clustered and the brown module was significantly correlated with cluster 3 (r=0.46, p=2e-06) (Figure 3B and C). Therefore, the brown module was chosen as the focus module for the research. The GO pathway showed most of the genes in this module to be related to TME, including protein polyubiquitination, interferon-gamma-mediated signalling pathway and autophagosome (Figure 4A). The KEGG pathways suggested ubiquitin-mediated proteolysis, bacterial invasion of epithelial cells, Kaposi sarcoma-associated herpes virus infection and Chemokine signalling pathway may be activated (Figure 4B). On the STRING, all genes in the brown module were then used for building a PPI network and an interaction score greater than 0.9 was taken as the threshold (Figure 5A). The degree scores for each gene were determined by the Network Analyzer in Cytoscape. A degree value ≥20 was defined as a hub gene and 30 hub genes were screened in total (Figure 5B). At the same time, the TCGA and ICGC databases were used for analysing the relationship between hub gene expression and HCC stage (Figure 6A and B). GRB2, SH3KBP1, LCK and UBE2D1 expressions were all found to be correlated with HCC stage in both databases. Furthermore, a high expression of GRB2 was found to have the ability to predict poor prognosis in advanced HCC patients (Figure 6C and D) and higher pathological stage (Figure 6E and F) in both databases. These results suggest that GRB2 expression has a close correlation with HCC patient stage and prognosis. Therefore, the role of GRB2 in HCC requires further analysis.

TME has been found to play a crucial role in immunotherapy in recent studies.23,24 Therefore, the relationship between GRB2 and immune cells was examined (Figure 7A and B). The results showed the high expression of GRB2 to be positively related to immune cells, including CD8+T cells, central memory CD8+T cells, activated CD4 + T cells, activated CD8 T cells and Myeloid dendritic cells. The single-cell sequencing also found there to be a significant increase in tumour suppressive immune cells, including CD4 Conv T cells, Treg cells, T profile cells, CD8 +T cells and CD8 Tex cells infiltration in the gene with a high expression of GRB2 from GSE140228_10X and GSE98638 (Figure 8A and C). The results suggest that GRB2 may affect the prognosis of advanced HCC patients by altering the immune cell infiltration level.

Somatic mutations may serve as immunotherapy response predictors,25 so 87 HCC patients were processed with gene mutation data using the maptools R package. The most common type of DNA base mutation in both low and high expression groups was found to be C>T (Figure 9A and B). The differences in somatic mutations between high and low expression groups were also compared. The most commonly mutated genes in the high expression group were TP53 (35%), TTN (26%), CTNNB1 (26%) and MUC16 (23%) and the most commonly mutated genes in the low expression group were CTNNB1 (27%), TTN (20%), APOB (16%) and TP53 (16%). Furthermore, the median mutation number in the low expression group was 70.5 and 64 in the high expression group (Figure 9C and D). Genes with different mutations may contribute to immune infiltration characteristics.

Several studies have found that TMB and immune infiltration can affect the treatment effect of HCC patients.26–28 Therefore, the relationship between GRB2 expression and HCC treatment was further evaluated. The TIDE algorithm was used for verifying the capacity of GRB2 expression in ICIs. The results found that most HCC patients who had a better response to ICIs had higher GRB2 expressions (p = 0.0065) and higher TIDE scores in comparison to a low expression of GRB2 (90%VS 80%) (Figure 10A). Conversely, patients with low GRB2 expressions had a higher response rate than those with a high GRB2 expression to sorafenib (42% vs 21%) and TACE (70% vs 41%) (Figure 10B and C). The results demonstrate that GRB2 has a high value in predicting ICI, sorafenib and TACE treatment efficacy among advanced HCC patients.

The potential biological functions of GRB2 in HCC were explored preliminarily using GSEA. From GO, it was found that a high GRB2 expression is often enriched in immune-related pathways, such as B cell proliferation, cytokine activity, cytokine binding and cytokine receptor binding (Figure 11A and B). KEGG suggested that GRB2 is involved in the interferon-related pathways, such as JAK-STAT signalling pathway and cytosolic DNA sensing pathway (Figure 11C and D). These results may further help with the elucidation of the potential mechanism of action of GRB2 in HCC.

Further confirming the clinical significance of GRB2 expression in HCC, the RNA expression level of GRB2 was found to be upregulated in 24 of the 30 pairs of HCC tissue and adjacent non-cancer tissues (80%) (p < 0.05) (Figure 12A). Similarly, IHC and Western blot were performed as a means of confirming that GRB2 was highly expressed in HCC protein levels (Figure 12B, C and F). In addition, immunohistochemical chip results showed there to be a positive correlation between GRB2 expression and tumour staging (Figure 12D and E). GRB2 expression was found to have a significant correlation with the presence or absence of liver cirrhosis (p = 0.007), tumour size (p = 0.002) and TNM stage (p = 0.049) (Table 1). The overall survival curve and disease-free survival curve can be seen in Figure 12G. The results of this study matched those above.

TME is an important factor that affects the efficacy of sorafenib, TACE and ICIs among advanced HCC patients.29,30 Recent studies have shown that ubiquitination-related enzymes play a crucial regulatory role in tumourigenesis and TME.31,32 Therefore, screening out ubiquitinated genes that regulate key molecules of immune infiltration among advanced HCC patients may provide new targets for improving HCC treatment efficacy. This study used ssGSEA for scoring immune infiltration in 90 advanced HCC patients from the TCGA database and divided them into three clusters using unsupervised cluster analysis. In comparison to the other clusters, cluster 3 was found to be the most immune enrichment type. WGCNA was then established as a means of obtaining the most relevant ubiquitinated gene module for immune infiltrating cells and PPI was performed for identifying hub genes in this module. The GRB2 were screened out, which was correlated with immune cell infiltration, prognosis and HCC stage in both the TCGA and ICGC databases. Growth factor receptor-bound protein 2 (GRB2) is a scaffold protein that has two SH2 and one SH3 domain. It has been proven that the SH3 domain regulates ubiquitination modifications of proteins by binding to ubiquitin.33 In addition, several studies have reported that GRB2 contributes to the progression or metastasis of several carcinomas, including colorectal, gastric and breast cancer.34–36 With liver cancer, lncRNA AC092171.4 upregulates GRB2 and promotes HCC progression by competitively binding miR-1271.4.37 Yu et al found the inhibition of GRB2 expression significantly inhibits cancer cell proliferation and survival.38 Furthermore, this work provides evidence that the mRNA and protein expression levels of GRB2 are high in HCC tissues and HCC patients with high expressions of GRB2 have worse prognosis. It was also found that the expression of GRB2 has a positive correlation with immune infiltration cells (CD4 Conv T cells, Treg cells, T profile cells, CD8+ T cells and CD8Tex cells), which plays a central role in HCC immunotherapy and improves patient prognosis.39,40 Recent studies have shown that PD1 or CTLA-4 inhibitors have obvious effects on Treg cells mediated immunosuppression and anti-tumour immunity.41 With NSCLC, patients with high expressions of PD1 and CD8+T have a longer progression free survival (PFS) after they receive ICIs than those with low expressions of PD1 or CD8+T.42 Furthermore, GRB2 plays a crucial role in immune infiltration cells. DNaM-1 promotes NK cell activation by coupling GRB2 and GRB2 regulated CD8+T cells and promoting the secretion of IL-2.43,44 Tumour mutation burden (TMB) also affects response to HCC immunotherapy. Mutation analysis has shown that the immune-related mutation rate is significantly higher in patients with high expressions of GRB2 expression than those with low expressions and the types of mutated genes differ widely between the groups. TP53 is the most frequently mutated gene in the high expression group (35%). With breast cancer, TP53 mutation significantly promotes the infiltration of Tregs, T helper cells and M0 macrophages.45 TP53 mutation has an inverse correlation with anti-CTLA-4 efficacy in melanoma.46 Therefore, TMB is an important prognostic factor for HCC and TMB-related genes could be potential HCC therapeutic targets. In the individual patient-level analysis of patients who were treated with ICIs, TERT, CTNNB1, BRD4, MLL or TP53 mutations were found to be associated with a significantly worse survival rate.47 Recent studies have also found that cabozantinib and nivolumab combination therapy improves the prognosis of advanced HCC patients with high TMB significantly.27 Sorafenib, TACE and ICIs are currently the main treatment modalities for advanced HCC patients.48 ICI efficacy is directly proportional to high GRB2 expressions. Conversely, patients with high GRB2 expressions have poor efficacy on sorafenib and TACE. Recent studies have reported that miR-378a-3p targets MAPK1/GRB2 sensitise ovarian cancer cells to cisplatin.49 Lymphocylin inhibits the GRB2/ERK/STAT3 signalling pathway and reduces the resistance of NSCLC cells to EGFR-TKI.50 Another study has reported that patients with high GRB2 expressions have a better response to PD-1/PD-L1 and TP53 frameshift mutation, which is potentially related to the high immune cell infiltration that is caused by GRB2, and this could improve the efficacy of PD-1/PD-L1.51 The above results provide further proof that GRB2 may exert critical function in HCC treatment. Finally, GSEA was performed as a means of exploring the mechanism of GRB2. The results suggest that GRB2 is involved in the cancer-related and cellular regulation pathways, such as the JAK-STAT3 signalling pathway, cytosolic DNA sensing pathway and regulation of autophagy. Recent studies have shown that the overexpression of STAT3 results in the up-regulation of Mcl-1 expression, which makes the hepatoma cell Huh7 resistant to sorafenib.52 With prostate cancer, JAK2-STAT3 inhibits the cGAS-STING and plays a crucial role in tumour rejection and immune cell infiltration.53 With liver cancer, a combination of sorafenib and STAT3 knockout induces endoplasmic reticulum stress and activates the cGAS-STING axis.54 GRB2 also regulates the genomic stability of (DNA damage repair) DDR by affecting the Rad51 expression.55 Therefore, it is supposed that GRB2 may activate STAT3 and the cGAS-STING axis through DDR and affect anti-tumour immunity. Recent studies have determined that a combination of STAT3 antiagonists and STING agonists could inhibit tumour growth among mice by increasing CD8+ T cells while reducing Tregs and MDSCs.53 Furthermore, STING agonist cannot individually induce IFN expression in prostate cancer cell DU145, but inhibition of the JAK2-STAT3 pathway may restore responsiveness to STING agonist.

To summarise, through a comprehensive evaluation of the molecular, cellular and clinical characteristics of GRB2, it was discovered that there may be important implications for the selection of sorafenib, TACE and ICI strategies for advanced HCC patients. However, there are limitations to the study that require refinement. Further experiments are required exploring the potential biological mechanism of GRB2 in HCC treatment.