TMEM106C, BSG, COPE, CDCA8, KPNA2, LIG1, UQCRH, and CCT5: Predictive of Survival and Immunotherapy Resistance in Hepatocellular Carcinoma

The following data were collected for our analyses: 1The scRNA‐seq data of HCC containing three HCC and three healthy liver samples were downloaded from the dataset GSE162616↗ [15].2TCGA was accessed to download the data of clinical information and gene expression of patients with liver hepatocellular carcinoma (LIHC). The FPKM value in the RNA‐seq data was transformed to transcripts with TPM and the log₂ transformation was carried out. Samples with survival > 30 days were retained, and 342 HCC samples and 50 control samples were accordingly collected.3The clinical information of HCC patients and the gene expression data were obtained from the dataset GSE43619↗ (platform: GPL10558↗). A total of 88 tumor samples were included [16].4The ICGC‐LIRI‐JP cohort (http://lifeome.net/database/hccdb/home.html↗) was downloaded from the HCCDB database [17]. A total of 211 HCC and 177 control samples were incorporated.

The “Read10X” function of the “Seurat” R package was applied to read the scRNA‐seq data of each sample [18, 19]. Cells with gene count between 200 and 3000 and with mitochondrial gene count < 10% were retained and standardized using the “SCTransform” function. Following the dimensionality reduction via principal component analysis (PCA), the intersample batch effects were removed using the “Harmony” R package [20]. The “RunUMAP” function was applied for dimensionality reduction, and the “FindNeighbors” and “FindClusters” functions (parameters: dims = 1 : 25 and resolution = 0.1) were employed for clustering the cells. The obtained cell clusters were finally annotated using the markers in the CellMarker2.0 database.

Additionally, the “AUCell” R package was applied to calculate the senescence score in different cell populations [21].

The data from TCGA were applied for the analysis to reveal the DEGs with the help of the “limma” R package [22]. The relevant DEGs were sorted at the parameters |log2FC| > 1 and adjusted p value < 0.05. The effects of DEGs were explored according to the three terms of gene ontology (GO) enrichment analysis.

The prognostically relevant genes were firstly sorted via univariate Cox regression analysis and reduced using LASSO Cox regression analysis by the “glmnet” R package [23]. The key genes and their coefficients were obtained via stepwise regression and the risk score for the model was calculated using the Formula (1). (1)Riskscore=Σβi×Expi

(β in the formula is the coefficient number, and i represents the gene level of expression.)

The patients were then stratified into high‐ or low‐risk groups based on the median of the calculated risk score, and the Kaplan–Meier (KM) curve was plotted for survival prognostication. The relevant receiver operator characteristic (ROC) curve was plotted, and the AUC value was calculated based on the “timeROC” R package [24].

The immune cell infiltration degree was estimated using two complementary algorithms, “MCPCounter” and “TIMER”, to ensure robust and comprehensive profiling of the tumor microenvironment. The scores of 28 types of tumor‐infiltrating lymphocytes (TILs) were further calculated using the “GSVA” R package [25], which enables a sensitive, nonparametric assessment of immune cell activity based on gene set variation analysis, thereby providing a detailed landscape of lymphocyte subsets within the samples.

The difference on the predicted immunotherapy response in patients of the two groups was compared using the “TIDE” algorithm. Meanwhile, the level of immune checkpoint–relevant genes in the two groups was also compared, and the corresponding results were visualized in a heatmap.

Human immortal adult liver epithelial cell line THLE‐2 (C5664, RRID: CVCL_3803) and HCC cell line HuH7 (C5176, RRID: CVCL_0336) were all purchased from BD Bio (Hangzhou, China) and cultured as follows: THLE‐2 cells were grown in bronchial epithelial cell growth medium (CC‐3170, Lonza, Basel, Switzerland) containing 10% fetal bovine serum (FBS, F814‐500, BD Bio, China), epidermal growth factor (5 ng/mL, 01‐107, Merck Sigma, Darmstadt, Germany), and Phosphoethanolamine (70 ng/mL, P0503, Merck Sigma, Germany). HCC cells HuH7 were grown in high‐glucose DMEM (L1004‐500, BD Bio, China) supplemented with 10% FBS. All cells were identified via short tandem repeat analysis and confirmed negative for mycoplasma contamination, which were incubated in an incubator at 37°C with 5% CO₂.

Based on the results and the existing studies [26], CDCA8 was applied as the target gene of interest for the knockdown assay. Accordingly, the small interfering RNA targeting CDCA8 and the corresponding negative control with scramble target sequence were all customized and ordered from GenePharma (Shanghai, China), which were then transfected into HCC cells HuH7 using Lipofectamine 2000 transfection reagent (11668027, Invitrogen, Carlsbad, CA, United States) as per the manuals. The target sequence of the small interfering RNAs was listed in Table S1.

The TriZol reagent (15596‐026, Invitrogen, United States) was employed to isolate total cellular RNA as per the manuals and the concentration of isolated RNA was thereafter quantified in a spectrophotometer (ND‐2000, ThermoFisher Scientific, Waltham, MA, United States). A commercial 1st strand cDNA synthesis kit (D7170S, Beyotime, Shanghai, China) was applied to synthesize the cDNA from 1 μg of total RNA, and the PCR was then performed using SYBR Green qPCR Mix assay kit (D7260, Beyotime, China) and CFX384 real‐time PCR detection system (Bio‐Rad, Hercules, CA, United States) with the indicated primers (detailed information is available in Table S2). The relative mRNA level was gauged with the 2^−ΔΔct method using GAPDH as internal control [27].

The transfected HCC cells HuH7 were seeded in a 96‐well plate at the concentration of 2 × 10³ cells per well and cultured for 24, 48, and 72 h, following which 10‐μL CCK‐8 cell viability test solution from the assay kit (C0037, Beyotime, China) was added for an additional 4‐h culture. Thereafter, the OD value was read using a microplate reader (iMark, Bio‐Rad, United States) at 450 nm.

The migration of HCC cells HuH7 following the transfection was explored via scratch assay. Concretely, the transfected cells were grown in a 6‐well plate at the concentration of 1 × 10⁵ cells/well and cultured overnight. Thereafter, a 200‐μL pipette tip was applied to make a scratch on the monolayer once cells were completely confluent, and cells were further incubated in an incubator at 37°C for 48 h. An inverted optical microscope (IM‐300, Optika Microscopes, Ponteranica, Italy) was adopted to observe the migrated cells and the degree of wound closure (%) was additionally evaluated.

Prethawed matrix gel (C0372, Beyotime, China) was added to the top transwell chamber (pore: 8 μm, code. 3422, Corning, Inc., Corning, NY, United States) beforehand and 2 × 10⁴ transfected HuH7 cells were populated in the chamber with 200‐μL serum‐depleted media, whereas the bottom transwell chamber was filled with 750‐μL complete culture media with 10% FBS. Following the culture period of 48 h, a cotton swab was applied to remove the cells in the upper chamber, and those cells invaded to the bottom chamber were serially fixed in 4% paraformaldehyde (P0099, Beyotime, China) for 10 min and stained in crystal violet (C0121, Beyotime, China) for another 10 min. All invaded cells were finally observed under an inverted optical microscope to quantify the number as needed.

All computational analysis was realized using R software (Version 3.6.0), and all data of laboratory assays were processed via GraphPad Prism software (Version 8.0.2). The Wilcoxon test was applied to compare the difference in two continuous variables, and the log–rank test was adopted to compare the survival of patients in different groups. For the experimental data, unpaired t‐test, one‐way, and two‐way ANOVA were applied to compare the data. All data with statistical significance were denoted by the asterisks at the threshold of p < 0.05, and those without statistical significance were marked with “ns” (nonsignificant, p > 0.05).

The data of scRNA‐seq were analyzed in the beginning to reveal some SRGs in HCC. Following the procedures of cell filtering, standardization, dimensionality reduction, and clustering, 80,997 cells were identified and allocated to eight main clusters with their distinct markers (Figure 1a,b). The calculation on the percentage of these eight cell clusters in HCC and healthy samples suggested that the percentage of NK cells was evidently higher in HCC samples (Figure 1c,d). Thereafter, all NK cells were extracted and reclustered into five main populations. The senescence score of these populations was additionally calculated, revealing a relatively higher score of Population 4 (Figure 1e,f). These findings revealed the specific cell populations for our research and paved the way for our additional analyses on the genes of interest in HCC.

The senescence score of samples in TCGA was calculated via ssGSEA, and a higher senescence score was seen in HCC samples than that in healthy samples (Figure 2a). Then, the HCC samples were allocated to the high or low score group by the median value of the senescence score, and the survival of patients in these two groups was compared. A poor prognosis was noticed in the high score group (Figure 2b), and the senescence score gradually increased with the severity of the clinical stages (Figure 2c,d). To validate our discoveries, similar analyses were also implemented on the ICGC‐LIHC‐JP cohorts, and similar results were seen as well. In other words, the senescence score was higher in HCC samples, and a higher senescence score was indicative of a worse prognosis (Figure 2e,f), thus demonstrating the involvement of senescence in HCC.

The DEGs of HCC and normal samples in TCGA were firstly analyzed and the relevant results were displayed in a Volcano plot (Figure 3a). The analyzed DEGs were then intersected with the genes highly expressed in Population 4 of NK cells, and 222 common genes were obtained (Figure 3b). GO enrichment analysis was then initiated on these common genes. The corresponding results have manifested that these genes were mainly enriched in DNA conformation change and nuclear division, chromosome segregation (Figure 3c). These findings thus hint at the potential modulatory effects of these overlapped genes on cell cycle and cell division so as to participate in the progression of HCC.

The univariate Cox regression analysis was then performed on the 222 common genes to sort the prognostically relevant genes. Further, LASSO regression analysis and stepwise regression were applied to reduce the number of these prognostically relevant genes. Finally, eight feature genes were identified and applied to construct the risk model using the formula below (Figures 4a,b):

Risk score = 0.301∗TMEM106C + 0.374∗BSG + (−0.5∗COPE) + 0.341∗CDCA8 + 0.343∗KPNA2 + (−0.776∗LIG1) + 0.282∗UQCRH + 0.307∗CCT5

The risk score of each sample in TCGA cohort was calculated, and the samples were allocated based on the median of the calculated risk score to the high‐ or low‐risk group. The efficacy of the risk model was then validated in the cohorts of TCGA‐LIHC. Based on the results, the poorer prognosis was observed in patients of high‐risk group than that of low‐risk group (Figure 4c). Additional validation using ICGC‐LIHC and GSE43619↗ cohorts has shown that high‐risk group of patients tend to have a shorter overall survival (Figure 4d,e). The ROC curve and the AUC values were further plotted and calculated. In TCGA‐LIHC, the AUC value of the risk model in predicting the 1‐, 2‐, 3‐, 4‐, and 5‐year survival was 0.87, 0.75, 0.72, 0.71, and 0.75 (average: 0.76, Figure 4f). As to the ICGC‐LIHC‐JP cohorts, the AUC value of the risk model in predicting the 1‐, 2‐, 3‐, and 4‐year survival was 0.65, 0.73, 0.75, and 0.74 (average > 0.7, Figure 4g). Besides, in the dataset GSE43619↗, the AUC value of the risk model in evaluating the 1‐, 2‐, 3‐, 4‐, and 5‐year survival was 0.89, 0.9, 0.78, 0.61, and 0.57 (average: 0.75, Figure 4h). These discoveries hence demonstrated the efficacy of risk score in stratification of HCC patients into different risks.

Furthermore, a series of laboratory assays were carried out to validate the potential involvement of the eight feature genes in HCC. Based on the data from quantification test in HCC cells HuH7 and human immortal adult liver epithelial cell line THLE‐2, the differential expression of the eight feature genes in HCC cells was observed, with a profound elevated expression of TMEM106C, CDCA8, KPNA2 and UQCRH in HCC cells (Figure S1a). Given that CDCA8 is significantly overexpressed in HCC cells and has been clearly documented in the literature to promote mitotic chromosome segregation and HCC progression [26], it was selected for subsequent functional validation experiments. We validated the efficiency following CDCA8 knockout (Figure S1b). Subsequently, according to the results of CCK‐8 (Figure S1c), scratch (Figure S1d), and Transwell assays (Figure S1e), the knockdown of CDCA8 via small interfering RNA could repress the viability, migration, and invasion of HuH7 cells in vitro (Figures S1c, S1d, and S1e).

To improve the prognostic assessment of risk score in clinical practice, both univariate and multivariate Cox regression analyses were carried out as needed. The results have manifested that both clinical stage and risk score were the independent prognosis factors (Figure 5a,b), hinting that these two factors may have a profound impact on the survival of HCC patients. Thereafter, a nomogram incorporating the clinical stage and risk score was established (Figure 5c). The corresponding calibration curve and decision curve analysis showed that the nomogram is close to the actual prediction on the 1‐, 3‐ and 5‐year survival (Figure 5d) and that the integration of clinical stage and risk score can provide more support on clinical decision than the application of clinical stage or other prognostic model alone (Figure 5e), thus demonstrating the potential of nomogram on precise prognostication in HCC.

The immune infiltration status of TCGA‐LIHC data was analyzed in the following three algorithms to reveal the potential effects of tumor microenvironment (TME) on the prognostication of HCC. Firstly, the immune infiltration score of each sample in TCGA‐LIHC was determined using ssGSEA, revealing that the score was evidently higher in the samples of high risk. In particular, the infiltration degree of myeloid‐derived suppressor cells (MDSC), T cells (regulatory T cells, activated CD4⁺ T cells, central memory CD8⁺ T cells), and dendritic cells (activated dendritic cells and immature dendritic cells) was evidently higher in the samples of high risk (Figure 6a).

Secondly, the infiltration status of TILs was evaluated using the TIMER algorithm, and the corresponding results have manifested that the score of TILs in the samples of high risk (Figure 6b). Additionally, MCPCounter analysis suggested the high infiltration of TILs in the samples of high risk (Figure 6c). These evidences collectively demonstrated the high immune infiltration in HCC patients of high risk.

The TIDE algorithm was employed to evaluate the therapeutic response of HCC patients to the immune checkpoint inhibitors (ICIs) using the data from TCGA‐LIHC. It was illustrated that the percentage of patients at high risk to ICIs was evidently lower than that of low risk (Figure 7a), suggesting the poor response to immunotherapy. Further analyses have unveiled that the TIDE score and the immune evasion‐related score (exclusion score) were higher in HCC patients at high risk (Figure 7b,c). Such discoveries hinted at the possible involvement of immune evasion in these patients at high risk. Additionally, comparison of the expression levels of immune checkpoint–relevant genes showed higher expression levels of these genes in HCC patients at high risk (Figure 7d).