Key Candidate Prognostic Biomarkers Correlated with Immune Infiltration in Hepatocellular Carcinoma.

BACKGROUND: Hepatocellular carcinoma (HCC) is the most common subtype of primary liver cancer, which causes ~800,000 deaths annually world-wide. Immune checkpoint inhibitor (ICI) has reformed cancer therapy and achieved unprecedented results in various malignancies, including HCC. However, the response rate of immunotherapy is very low in HCC. Considereing the complicated and unique immune status in liver, we hypothesize that critical molecules will affect prognosis and correlate with immune context in the tumor microenvironment of HCC.
METHODS: Using Kaplan-Meier plotter, GEPIA2 and Integrative Molecular Database of Hepatocellular Carcinoma (HCCDB), survival genes and their prognostic value were estimated in HCC. Based on Tumor Immune Estimation Resource (TIMER), association between survival genes and immune infiltration was examined in HCC. FunRich and STRING were used to analyze gene ontology and protein-protein interaction (PPI) Network, qRT-PCR was used to measure mRNA level of candidates; and a Cell Counting Kit-8 was used to measure proliferation of HCC cell line.
RESULTS: Using multiple databases, we identified 36 key prognostic genes highly expressed in HCC and associated with poor survival of patients. Meanwhile, the 36 gene signatures correlated with immune infiltration in HCC. Moreover, these genes were significantly associated with exhausted T cells and polymorphonuclear myeloid-derived suppressor cells (PMN-MDSCs) in HCC. Among the 36 key genes, SKA3, SGOL2, SPINDOC, TEDC2, TMCO3 and NUP205 were highly expressed in tumor samples compared with adjacent normal tissues in our HCC cohort (n=22). Additionally, proliferation of SMMC7721 cell line was inhibited when it interfered with SiRNA of each gene.
CONCLUSION: The 36 genes may serve as potential prognostic biomarkers and molecular targets to ameliorate tumor immune microenvironment (TIME) in HCC and therefore represent a novel avenue for individualized immunotherapy in HCC.

Introduction

Hepatocellular carcinoma (HCC), which accounts for approximately 90% of primary liver cancers, is the sixth most common malignancy and the second leading cause of cancer-related mortality in the world.1,2 Although treatments with Sorafenib and Regorafenib show a modest survival benefit, their overall anti-tumor efficacies are still limited.3–6 Over recent decades, the landscape of cancer treatment has been revolutionized by immunotherapies.7 While such therapies as checkpoint blockade can lead to unprecedented clinical benefits for malignant melanoma, non-small cell lung cancer (NSCLC), bladder cancer and Hodgkin lymphoma, etc.,7–11 the response rate is low, and it is much lower in HCC,12 the underlying mechanisms remain to be elucidated. As a major immune organ, liver acts as a central immune-modulator to protect other organs and systems by maintaining immunologic tolerance.13 Sophisticated and unique immune sets in liver and HCC tumor microenvironment (TME) play important roles in virus infection, inflammatory reaction, immune evasion and response to immune therapy.13 Recently, various studies have identified that tumor-infiltrating lymphocytes in TME are associated with prognosis and immunotherapy response of HCC patients. For example, increased regulatory T cells (Treg) correlate with poor overall survival and resurgence-free time, while increased ratio values of B/T cells in TME correlate with longer survival in HCC.14,15 In addition, immune cell infiltration of the primary tumor affects individual response to checkpoint inhibitors in metastatic melanoma independent of PD-L1 status.16 Thus, there is an urgent need to uncover the molecular mechanisms of immune infiltrates in TME and to identify novel immune-related therapeutic targets in order to improve immunotherapy in HCC.

With rapid advances in next-generation sequencing (NGS) and bioinformatics, accumulated data have helped us to understand complicated biological characteristics of tumors in multiple dimensions. Simultaneously, more and more databases have been generated to analyze cancer pathogenesis comprehensively, such as Gene Expression Omnibus (GEO, ) and The Cancer Genome Atlas (TCGA, ),17 which have already molecularly characterized over 20,000 primary tumors and their matched normal samples among 33 types of cancer. According to the abundant data from TCGA and other databases, various related online visualized analyzing tools and software have also been created, such as GEPIA2 (), which focuses on RNA sequencing data of 9736 tumors and 8587 normal samples from TCGA and GTEx projects. This can provide customizable functions, including differentially expressed survival-related genes in various cancers.18 Another tool, called Tumor Immune Estimation Resource (TIMER, ), is a comprehensive resource for systematic analysis of immune infiltrates across multiple cancer types.19,20 Such tools based on comprehensive databases extensively help us to uncover new potential prognostic markers in various cancer types as well as complicated associations between molecules and tumor immune microenvironment (TIME), so as to provide new potential targets for HCC immunotherapy, and to improve the precise treatment in HCC patients.

In our present study, we identified 36 key prognostic genes highly expressed in HCC and associated with poor survival of patients through multiple databases. Meanwhile, the 36 gene-signature correlated with immune infiltration in HCC. Moreover, these genes were significantly associated with exhausted T cells and polymorphonuclear myeloid-derived suppressor cells (PMN-MDSCs) in HCC. Among them, SKA3, SGOL2, SPINDOC, TEDC2, TMCO3 and NUP205 were highly expressed in tumor samples compared with adjacent normal tissues in our HCC cohort (n=22), Additionally, proliferation of SMMC7721 cell line was inhibited when it interfered with SiRNA of each gene.

Materials and Methods

Gene Analysis in GEPIA2

GEPIA2 is an updated version of GEPIA that analyzes RNA sequencing expression data of 9736 tumors and 8587 normal samples from TCGA and GTEx projects, using a standard processing pipeline. It provides customizable functions such as tumor/normal differential expression profiling, according to cancer types or pathological stages, patient survival analysis, similar gene detection, correlation analysis and dimensionality reduction analysis ().18 In our study, we applied survival analysis to the most differential survival genes in hepatocellular carcinoma (HCC). We analyzed correlation between these genes and immune infiltrates in HCC through TIMER. We defined overall survival as prognosis and the median as group cutoff.

Correlation Analysis Through TIMER Database

Tumor Immune Estimation Resource (TIMER) is a comprehensive resource for systematic analysis of immune infiltrates across diverse cancer types ().19,20 Relative abundances of 6 immune infiltrates (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages and dendritic cells) are estimated, which will be validated using pathological estimations.19 In our present study, we used TIMER to analyze association between the most differential survival genes and immune infiltrates in HCC. Moreover, associations between each gene and various markers of tumor-infiltrating immune cells were analyzed, including CD8+T cells, T cells (general), B cells, monocytes, TAMs, M1 macrophages, M2 macrophages, neutrophils, natural killer (NK) cells, dendritic cells (DCs), T-helper 1 (Th1) cells, T-helper 2 (Th2) cells, follicular helper T (Tfh) cells, T-helper 17 (Th17) cells, Treg, exhausted T cells, myeloid-derived suppressor cells (MDSC) and mucosal-associated invariant T (MAIT) cells, as referenced in prior studies.13,41,46,51–56 Then, Spearman correlation and estimated statistical significance of a pair of defined genes in HCC were generated through correlation modules. The targeted gene was applied to x-axis, and related markers were represented on y-axis. Relatively, gene expression level was displayed with log2 RSEM.

Gene Expression and Prognostic Analysis Through HCCDB Database

Integrative Molecular Database of Hepatocellular Carcinoma (HCCDB) was contributed by Tsinghua University and National Center for Liver Cancer & Shanghai Eastern Hepatobiliary Surgery Hospital. This database curated 15 public HCC expression datasets with up to ~4,000 clinical HCC samples and serves as a one-stop online resource for HCC gene expression and other comprehensive analyses, including prognostic and co-expression analysis ().21 In our study, the most differential survival genes associated with immune infiltrates were analyzed by HCCDB, to compare the expression levels of these genes in HCC tumors with normal tissues, as well as their prognostic value in HCC.

Kaplan–Meier Plotter Database Analysis

The Kaplan–Meier plotter is capable of assessing potential effects of 54,000 genes on survival in 21 cancer types. Primary purpose is meta-analysis based discovery and validation of survival biomarkers ().57,58 In our study, we used this database to confirm prognostic value of selected key genes in HCC. The HR with 95% CI was marked.

Gene Ontology and Protein–Protein Interaction (PPI) Network Analysis

For better understanding of key prognostic genes associated with immune infiltrates in HCC, cellular component, molecular function and biological process were analyzed by Gene Ontology through FunRich software (). PPI network was analyzed by STRING online tool ().

SPINDOC, TEDC2, NUP205, SGOL2, SKA3 and TMCO3 Expression Analysis

This study was approved by Ethical Committee of Beijing Shijitan Hospital and performed by complying with all relevant ethical regulations [No. sjtkyll-lx-2021(22)]. Written informed consent from each patient was obtained prior to study commencement. All processes were conducted in accordance with the Declaration of Helsinki.

About 10 mg HCC sample or adjacent normal tissue was disrupted and homogenized in 1 mL RNA-Solv Reagent (OMEGA R6934) containing 2% 2-mercaptoethanol. Total RNA was isolated using the E.Z.N.A Total RNA Kit II (OMEGA R6934-02) according to the tissue protocol. ReverTra Ace qPCR RT Kit (TOYOBO FSQ-101) was used to synthesize the first-strand cDNA. Real-time qPCR was done on Quant Studio 6 Real Time PCR System with SYBR (Thermo Fisher Scientific). 18S was used to normalize relative expression of target genes.

About 1×10^5 SMMC7721 cells were seeded into 12-well plates. Lipofectamine 3000 reagent was added to transfect SiRNA (final concentration=0.1MM) into cells at 70–80% confluence, according to standard protocol (Thermo Fisher Scientific). Cells were incubated in incubator with 5% CO2 at 37°C. Two days later, total RNA was collected using E.Z.N.A Total RNA Kit II (OMEGA R6934-02) according to cell protocol, followed by qPCR, as described previously.

Proliferation Assays

Two thousand SMMC7721 cells were seeded in 96-well plates. One day later, SiRNA (final concentration=0.1MM) was transfected into cells as described above. Cell Counting Kit-8 (YESEN 40203, Shanghai, China) was used to measure proliferation from the 1st to 5th day after interfered with SiRNA. Cell viability (%) = [OD450 (SiRNA)-OD450 (BLANK)]/[OD450 (SiRNA-NC)-OD450 (BLANK)].

Statistical Analysis

The survival curves were generated by GEPIA2, Kaplan–Meier plotter and HCCDB. Spearman correlation was used to evaluate association between gene expression and a specific parameter. The strength of correlation was defined by absolute value as follows: 0.00–0.19 “very weak”, 0.20–0.39 “weak”, 0.40–0.59 “moderate”, 0.60–0.79 “strong”, and 0.80–1.0 “very strong”. Statistical comparisons of mRNA levels between HCC tumors and adjacent normal tissues were conducted using paired t-tests or a one-way analysis of variance (ANOVA). Unpaired t-test was used to analyze viability of cells between NC and SiRNA group. A P-value of <0.05 was considered statistically significant.

Results

Identification of the Most Differential Survival-Related Genes and Their Association with Immune Infiltrates in HCC

Using TCGA database through GETIA2, the most differential survival-related genes were explored, which may serve as potential prognostic markers in HCC (Figure 1 and ). Five hundred genes were identified, including the top 10 genes HILPDA, CCDC58, B3GAT3, CLEC3B, CTB-147N14.6, SOCS2, LPCAT1, PES1, DNASE1L3 and AHSA1. Then, the association between each differential survival gene and immune infiltrates in HCC was analyzed through TIMER. We defined a relation between a gene and immune infiltration as existing when >4 types of immune cell (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages and dendritic cells) were associated with a target gene. Consequently, 119 of 500 genes were identified as being associated with immune infiltrates in HCC (Figure 1 and ).

Figure 1

Thirty-six key genes were identified from multiple databases. Differential genes: The 500 most differential survival genes in HCC through TCGA database; Immune infiltration: the differential survival genes which are generally correlated with immune infiltration in HCC in TCGA through TIMER; Key genes: the differential expression genes with prognostic value which are correlated with immune infiltration in HCC from the 500 most differential survival genes.

Identification of Key Candidate Prognostic Biomarkers Correlated with Immune Infiltrates in HCC

To identify new potential prognostic biomarkers in HCC, first we compared the expression levels of 119 genes () between HCC tumor and adjacent normal tissues, then we evaluated the prognostic values of these genes through the database Integrative Molecular Database of Hepatocellular Carcinoma (HCCDB), which has integrated 15 public datasets with up to 4000 clinical HCC samples and serves as a one-stop online resource for HCC gene expression and comprehensive analysis, including prognostic and co-expression analyses ().21 From the above 119 genes correlated with immune infiltrates, 36 genes with significantly different expression between HCC and adjacent tissues were identified (Figure 2 and ). Further analysis indicated that the different expression levels of the 36 genes between HCC and adjacent tissues exhibited notable prognostic value in HCC; a higher level of gene expression was associated with much poorer survival (Table 1). In addition, we also confirmed each gene’s prognostic value in HCC through Kaplan–Meier plotter database analysis () (Figure 3 and ). The results suggested that the 36 genes may serve as new potential prognostic markers in HCC.

Figure 2

The key prognostic genes were correlated with immune infiltrates in HCC. “partial.cor” means tumor purity-corrected partial Spearman correlation. (A) C11ORF84 (SPINDOC), (B) C16ORF59 (TEDC2), (C) CAPG, (D) CBX1, (E) CCNB1, (F) CDC20, (G) CDC45, (H) CDCA5.The immune infiltrates correlation analysis results of the other 28 genes can be found in with a highly yellow marker.

Table 1

Key Candidate Prognostic Biomarkers Correlated with Immune Infiltrates in HCC

Genes	Expression	Fold Change	Prognostic Performance
NUP205	UP	1.92	Unfavorable
CAPG	UP	7.13	Unfavorable
CBX1	UP	2.53	Unfavorable
CCNB1	UP	16.44	Unfavorable
CDC20	UP	30.06	Unfavorable
CDC45	UP	9.59	Unfavorable
CDCA5	UP	14.41	Unfavorable
CDCA8	UP	9.28	Unfavorable
CDK1	UP	19.15	Unfavorable
CDT1	UP	13.63	Unfavorable
CEP55	UP	12.43	Unfavorable
CLIC1	UP	3.32	Unfavorable
DLGAP5	UP	1.5	Unfavorable
E2F8	UP	25.33	Unfavorable
EXO1	UP	18.1	Unfavorable
EZH2	UP	2.85	Unfavorable
FABP5	UP	5.78	Unfavorable
G6PD	UP	3.44	Unfavorable
KIF11	UP	6.5	Unfavorable
KIFC1	UP	20.8	Unfavorable
MCM3	UP	3.77	Unfavorable
MCM7	UP	2.71	Unfavorable
MCM8	UP	2.45	Unfavorable
MKI67	UP	17.06	Unfavorable
NEK2	UP	29.18	Unfavorable
NUF2	UP	19.08	Unfavorable
OIPS	UP	1.21	Unfavorable
PRC1	UP	9.54	Unfavorable
SPINDOC	UP	2.4	Unfavorable
SGO1	UP	7.71	Unfavorable
SGO2	UP	6.92	Unfavorable
SKA3	UP	22.86	Unfavorable
SPC25	UP	13.19	Unfavorable
TEDC2	UP	5.92	Unfavorable
TMCO3	UP	2.88	Unfavorable
TRIP13	UP	7.38	Unfavorable

Notes: Up: The expression level of genes was higher in HCC compared with adjacent tissues in TCGA/HCCDB; Fold change: The expression level of genes in HCC compared with adjacent normal tissues in TCGA; Unfavorable: The higher expression level is correlated with poorer prognostic survival in TCGA/HCCDB.

Figure 3

Kaplan–Meier survival curves were used to confirm the prognostic value of the key genes from Kaplan–Meier plotter database. Overall survival (OS) was used in all the genes to evaluate their prognostic value. (A) C11ORF84, (B) C16ORF59, (C) CAPG, (D) CBX1, (E) CCNB1, (F) CDC20, (G) CDC45, (H) CDCA5.The OS prognostic value analysis results of the other 28 genes can be found in .

Correlation Between Key Biomarker Expressions and Immune Marker Sets in HCC

To our surprise, all 36 candidates were upregulated in HCC. We then applied TIMER to further explore molecular mechanisms underlying relationship between 36 key potential prognostic markers and TIME (Table 2, ). Interestingly, these 36 genes correlated with infiltration of CD8+T cells, general T cells, B cells, monocytes, tumor-associated macrophages (TAMs), and type 1 T helper cells (Th1). The strength of correlation varied from “very weak” to “moderate”. These genes showed no significant correlation with subtype of macrophages (M1/M2), natural killer cells, follicular helper T cells (Tfh), Th17 and mucosal-associated invariant T cells (MAIT). For neutrophils and dendritic cells (DCs), these genes showed significant correlation with ITGAM and ITGAX, respectively. For exhausted T cells, these genes were significantly positively correlated with PDCD1, CTLA4, HAVCR2 and TIGIT. These genes also showed general correlation with LAG3, BTLA, ENTPD1, CD27, TNFRSF9, WARS and CXCL13. Thus, these genes might be highly connected with T cell exhaustion in TME of HCC, which may therefore affect the survival of HCC patients.

Table 2

Correlation Between Each Key Gene and Immune Cell Marker Sets in HCC

		C11ORF84		C16ORF59		CAPG		CBX1		CCNB1
		Partial.Cor	Partial.P	Partial.Cor	Partial.P	Partial.Cor	Partial.P	Partial.Cor	Partial.P	Partial.Cor	Partial.P
CD8+T	CD8A	0.135	*	0.157	**	0.227	***	0.065	0.231	0.165	**
	CD8B	0.169	**	0.236	****	0.258	****	0.051	0.342	0.217	***
T (General)	CD3D	0.263	****	0.321	****	0.423	****	0.107	*	0.321	****
	CD3E	0.201	**	0.207	**	0.351	****	0.063	0.243	0.203	**
	CD2	0.204	**	0.222	***	0.405	****	0.052	0.337	0.230	***
B cell	CD19	0.241	****	0.244	****	0.294	****	0.226	***	0.224	***
	CD79A	0.221	***	0.167	**	0.270	****	0.182	**	0.142	**
Monocyte	CD86	0.294	****	0.254	****	0.491	****	0.193	**	0.285	****
	CSF1R	0.154	**	0.113	*	0.326	****	0.080	0.137	0.086	0.109
TAM	CCL2	0.146	**	−0.008	0.887	0.268	****	0.053	0.329	−0.089	0.098
	CD68	0.070	0.192	0.126	*	0.402	****	−0.022	0.678	0.183	**
	IL10	0.231	***	0.122	*	0.313	****	0.176	**	0.208	***
M1	NOS2	−0.125	*	−0.178	**	−0.098	0.070	−0.065	0.226	−0.198	**
	IRF5	0.133	*	0.168	**	0.253	****	0.179	**	0.170	**
	PTGS2	0.216	***	−0.055	0.311	0.247	****	0.116	*	−0.025	0.644
M2	CD163	−0.031	0.570	0.019	0.729	0.135	*	0.017	0.754	−0.022	0.681
	VSIG4	0.075	0.167	0.082	0.130	0.211	***	0.032	0.556	0.019	0.720
	MS4A4A	0.003	0.959	0.061	0.256	0.236	****	−0.053	0.325	0.021	0.694
Neutrophil	CEACAM8	0.092	*	0.060	0.112	0.008	0.838	0.065	0.088	0.129	**
	ITGAM	0.243	****	0.316	****	0.351	****	0.251	****	0.249	****
	CCR7	0.131	*	0.019	0.730	0.204	**	0.045	0.404	−0.012	0.830
Natural killer cell	KIR2DL1	−0.102	0.057	−0.091	0.091	−0.107	*	−0.108	*	−0.117	*
	KIR2DL3	0.071	0.190	0.063	0.244	0.026	0.637	0.109	*	0.058	0.284
	KIR2DL4	0.090	0.095	0.168	**	0.084	0.118	0.033	0.544	0.140	**
	KIR3DL1	−0.074	0.169	−0.078	0.147	−0.044	0.411	−0.021	0.694	−0.107	*
	KIR3DL2	0.054	0.313	0.059	0.275	0.032	0.550	0.087	0.106	0.011	0.845
	KIR3DL3	0.043	0.428	0.026	0.629	0.027	0.623	0.029	0.592	0.042	0.437
	KIR2DS4	−0.001	0.989	−0.036	0.500	0.002	0.967	−0.023	0.667	−0.040	0.464
DCs	HLA-DPB	0.119	*	0.146	**	0.337	****	−0.006	0.906	0.128	*
	HLA-DQB	0.097	0.071	0.164	**	0.278	****	0.009	0.870	0.139	**
	HLA-DRA	0.082	0.130	0.119	*	0.288	****	0.048	0.371	0.133	*
	HLA-DPA	0.108	*	0.089	0.099	0.285	****	0.024	0.654	0.107	*
	CD1C	0.140	**	−0.015	0.776	0.224	***	0.124	*	0.053	0.330
	NRP1	0.122	*	−0.251	****	0.069	0.204	0.193	**	−0.111	*
	ITGAX	0.329	****	0.268	****	0.551	****	0.180	**	0.289	****
Th1	TBX21	0.027	0.611	0.025	0.649	0.134	*	−0.044	0.417	0.000	0.995
	STAT1	0.166	**	0.091	0.092	0.284	****	0.234	***	0.216	***
	STAT4	0.190	**	0.172	**	0.349	****	0.081	0.133	0.161	**
	IFNG	0.188	**	0.255	****	0.238	****	0.145	**	0.282	****
	TNF	0.280	****	0.219	***	0.336	****	0.168	**	0.246	****
Th2	GATA3	0.227	***	0.130	*	0.325	****	0.088	0.105	0.149	**
	STAT5A	0.137	*	0.143	**	0.281	****	0.005	0.923	−0.012	0.827
	STAT6	−0.293	****	−0.301	****	−0.101	0.062	−0.245	****	−0.344	****
	IL13	0.083	0.123	0.075	0.164	0.041	0.446	0.043	0.422	0.072	0.183
Tfh	BCL6	−0.027	0.619	−0.187	**	−0.196	**	0.030	0.583	−0.176	**
	IL21	0.096	0.075	0.108	*	0.087	0.105	0.135	*	0.167	**
Th17	STAT3	−0.083	0.125	−0.225	***	−0.122	*	−0.043	0.427	−0.221	***
	IL17A	−0.029	0.591	−0.040	0.457	0.045	0.409	0.070	0.192	0.061	0.258
Treg	FOXP3	0.057	0.294	0.030	0.583	0.049	0.366	0.116	*	0.061	0.255
	CCR8	0.300	****	0.163	**	0.317	****	0.369	****	0.289	****
	STAT5B	−0.167	**	−0.323	****	−0.308	****	0.105	0.050	−0.349	****
	TGFB1	0.434	****	0.167	**	0.486	****	0.244	****	0.227	***
Exhausted T cell	PDCD1	0.329	****	0.286	****	0.365	****	0.188	**	0.340	****
	CTLA4	0.333	****	0.366	****	0.401	****	0.201	**	0.384	****
	LAG3	0.179	**	0.251	****	0.251	****	0.057	0.289	0.265	****
	HAVCR2	0.323	****	0.313	****	0.535	****	0.183	**	0.328	****
	GZMB	0.065	0.227	0.088	0.103	0.117	*	−0.055	0.307	0.071	0.187
	TIGIT	0.272	****	0.277	****	0.307	****	0.207	**	0.295	****
	BTLA	0.119	*	0.116	*	0.214	***	0.098	0.069	0.128	*
	CD244	0.037	0.495	0.075	0.165	0.118	*	−0.052	0.339	0.085	0.116
	CD160	0.064	0.236	0.140	**	0.044	0.419	−0.065	0.229	0.093	0.084
	ENTPD1	0.387	****	0.097	0.072	0.367	****	0.231	***	0.177	**
	C10ORF54	0.151	**	0.116	*	0.357	****	0.025	0.643	0.012	0.831
	CD27	0.252	****	0.295	****	0.345	****	0.138	*	0.263	****
	TNFRSF9	0.299	****	0.221	***	0.219	***	0.323	****	0.256	****
	MYO7A	−0.115	*	−0.109	*	0.033	0.539	−0.078	0.150	−0.130	*
	WARS	0.186	**	0.177	**	0.283	****	0.041	0.446	0.213	***
	CXCL13	0.197	**	0.230	***	0.193	**	0.148	**	0.240	****
	LAYN	0.142	**	−0.017	0.760	0.233	***	0.068	0.205	−0.006	0.914
	PHLDA1	−0.078	0.148	−0.323	****	−0.033	0.547	−0.045	0.400	−0.228	***
	SNAP47	−0.112	*	−0.001	0.980	−0.006	0.908	0.065	0.230	0.002	0.971
MDSC	CD14	−0.391	****	−0.261	****	−0.181	**	−0.466	****	−0.401	****
	FUT4	0.344	****	0.136	*	0.419	****	0.301	****	0.210	***
	CD33	0.169	**	0.139	**	0.285	****	0.121	*	0.161	**
MAIT	KLRB1	−0.001	0.985	0.025	0.644	0.208	***	−0.140	**	−0.047	0.387

Notes: * p<0.05, **p<0.01, ***p<0.001, ****p<0.00001.

Abbreviations: TAM, tumor-associated macrophage; Th, T helper cell; Tfh, follicular helper T cell; Treg, regulatory T cell; MDSC, myeloid-derived suppressor cells; MAIT, mucosal-associated invariant T cell; partial Cor, R value of Spearman correlation adjusted by purity; partial p, p-value of Spearman correlation adjusted by purity.

In addition, for Treg cells, these genes showed very weak correlation with FOXP3, were negatively correlated with STAT5B and positively correlated with CCR8 and TGF-B1, suggesting these genes could probably influence Treg differentiation in HCC TME. Moreover, for MDSCs, these genes were significantly and positively associated with FUT4 (CD15, the molecular marker of PMN-MDSC) but negatively correlated with CD14 (the molecular marker of M-MDSC). Thus, considering their important role in the TIME, the other underlying mechanism behind the prognostic roles of these genes in HCC may be due to their positive correlations with Treg and MDSC.

Key Biomarkers' Gene Ontology and Protein–Protein Interaction (PPI) Network Analysis

To better understand potential roles of 36 key prognostic genes associated with immune infiltration, FunRich software was used to analyze cellular component, molecular function and biological process. Most of these genes can be localized in nucleus (82.4%) and nucleoplasm (58.8%) (Figure 4A). The most enriched biological processes were protein binding (91.7%), chromatin binding (22.2%) and DNA replication origin binding (11.1%) (Figure 4B). Moreover, these genes were associated with various functions, including cell division (45.7%), chromosome segregation (20%), and DNA replication (20%) (Figure 4C). In addition, based on PPI analysis, there were significant interactions among 36 proteins (avg. local clustering coefficient=0.774, PPI enrichment p-value <1.0e-16, Figure 5).

Figure 4

(A) Functional enrichment analysis of key prognostic markers. (A) Cellular components, (B) Biological processes, (C) Molecular functions in gene ontology analysis using FunRich software. p<0.05 was considered a statistically significant difference.

Figure 5

Protein–protein interaction (PPI) network for key potential prognostic genes in HCC. Those 36 key prognostic genes related to the PPI network are analyzed by STRING online tool (avg. local clustering coefficient=0.774, PPI enrichment p-value<1.0e-16).

Expression of SPINDOC, TEDC2, NUP205, SGOL2, SKA3 and TMCO3 in HCC Samples

To confirm differential expression of key genes, we analyzed mRNA levels of SPINDOC, TEDC2, NUP205, SGOL2, SKA3 and TMCO3 in clinical samples. All these genes are highly expressed in HCC tumor tissues compared with adjacent normal tissues (Figure 6). The clinical information of the patients is summarized in .

Figure 6

Twenty-two paired HCC tissues and adjacent normal tissues were collected from Aug. 2011 to Dec. 2019 in our hospital. The qPCR results showed that all these 6 genes were more highly expressed in HCC samples compared with adjacent normal tissues.*p<0.05, **p<0.01, ***p<0.001.

Interfering mRNA Expression of SPINDOC, TEDC2, NUP205, SGOL2, SKA3 and TMCO3 Inhibit the Proliferation of the SMMC7721 HCC Cell Line

Considering that SPINDOC, TEDC2, NUP205, SGOL2, SKA3 and TMCO3 were highly expressed in HCC and higher expression is associated with poorer prognosis, we used specific SiRNA to knock down each gene respectively in human SMMC7721 HCC cell line. Interestingly, the proliferation of the SMMC7721 cell line was significantly inhibited when the expressions of SPINDOC, TEDC2, NUP205, SGOL2, SKA3 or TMCO3, respectively, were interfered with (Figure 7, and ).

Figure 7

CCK8 reagents were used to assay the proliferation ability after SiRNA interference. The most effective SiRNA from was used to interfere with the mRNA expression. *p<0.05, **p<0.01, ***p<0.001, ****p<0.0001.

Discussion

As the second leading cause of cancer-related death and the fifth most common cancer, primary liver cancer causes ~800,000 deaths each year, globally22. Generally, 80–90% of liver cancer cases are hepatocellular carcinoma (HCC).23 With advances in basic and clinical science, underlying mechanisms of HCC tumorigenesis and progression have been explored. The main risk factors for HCC are well known, including hepatitis B and C virus infection, alcohol consumption and aflatoxin.23 In addition, molecular pathogenesis of HCC was also revealed. Representative multiple targets drugs (Sorafenib and Regorafenib) have been widely applied in clinical practice and exhibit only a minor survival benefit, thus the prognosis of HCC remains very poor.3–6 Recently, immune checkpoint inhibitor (ICI) has reformed the paradigm of cancer therapy and achieved marvelous results in controlling several types of cancer. However, most cancer patients, especially those with HCC, show no response to ICI.12 Due to the complicated immune status in liver by itself, we hypothesize that some key molecules may affect its prognosis and correlate with the immune context in HCC.

In this study, we first analyzed the most differential survival-related genes in HCC using TCGA database and, through GEPIA2, 500 genes were identified (). We then analyzed the relationships between these genes and immune infiltrates in HCC through TIMER, and 119 genes were found to correlate with the immune infiltrates in HCC (Figure 1 and ). To further estimate their prognostic values in HCC, we first compared the expression levels of the 119 genes between HCC tumor tissues and adjacent normal tissues, and then evaluated their prognostic values in HCC through another database, HCCDB, which specifically collected data from HCC-related studies and comprehensively analyzed the characteristics of HCC.21 Among the 119 genes, 36 key genes exhibited significantly differential expression in HCC tumors compared with adjacent tissues. Meanwhile, the 36 genes showed unfavorable prognostic value in HCC. A higher expression of 36 genes correlated with much poorer survival (Table 1). In addition, the gene ontology analysis revealed that these 36 molecules were mostly located in the important organs of the cell (nucleus and nucleoplasm), participating in significant biological processes (protein binding, chromatin binding and DNA replication origin binding), and were associated with various key functions (cell division, chromosome segregation, DNA replication) (Figure 4). Moreover, (PPI network analysis indicated that these 36 molecules had mutual correlations (Figure 5). As we all know, uncontrolled cell division is one of the most significant characteristics of cancer and the results of our gene ontology and PPI analysis may partly confirm the significant correlation between these 36 genes and the prognosis of HCC.

Most of these genes (30/36) have been reported to correlate with the prognosis of HCC by other authors24 (), and our results were very much in accordance with those of previous studies. In addition, our results indicated that all 36 genes had a significant correlation with the immune infiltration in HCC. This kind of relationship, to the best of our knowledge, hads never been reported before. Furthermore, for the 6 critical potential prognostic genes (NUP205, SKA3, SPINDOC, SGOL2, TMCO3 and TEDC2) no type of relationship with HCC, including with the prognosis of HCC, to our knowledge, has ever been reported before. Moreover, the 6 genes were confirmed to highly express in HCC tumors in our own cohort (n=22) (Figure 6). Notably, the proliferation of SMMC-7721 human hepatocellular carcinoma cell line was inhibited when the expressions of NUP205, SKA3, SPINDOC, SGOL2, TMCO3 and TEDC2, respectively, were interfered with (Figure 7). Their underlying molecular mechanism is worth further exploration to discover potential therapeutic targets for HCC.

In a previous study, SPINDOC was shown to directly bind Spindlin1 (SPIN1) and to strongly disrupt the histone methylation reading ability of SPIN1.25 SPIN1 is a transcriptional coactivator with critical functions in embryonic development and plays emerging roles in cancer. However, the interaction between SPINDOC with SPIN1 in HCC is worth deep exploration.

The functions of TEDC2 remain unexplored. Using a triple-level (gene-, transcript-, and exon-level) analysis in lung adenocarcinoma transcriptomes with 77 paired tumor–normal tissues, one previous study proposed TEDC2 as a differentially expressed gene.26 So far, the role of TEDC2 in HCC is unclear, and fundamental and clinical investigation is needed.

NUP205 (nucleoporin 205) encodes a nucleoporin, which can actively transport proteins, RNAs and ribonucleoprotein particles between nucleus and cytoplasm. One previous study identified NUP205 as a TMEM209-interacting protein. Mass spectrometric analysis indicated that, when stabilizing NUP205, c-Myc was upregulated in the nucleus. Overexpression of TMEM209 and TMEM209–NUP205 interaction could drive the proliferation of lung cancer.27 It is promising to explore its role and underlying mechanisms in HCC tumorigenesis.

SGO2 (shugoshin 2), also known as SGOL2 or TRIPIN, plays an important role in cell cycle progression during meiosis.28,29 SGO2 contains a MAD1/CDC20-like–MAD2 interaction motif and competes with MAD1 and CDC20 for binding to Mad2.29 CDC20, one key prognostic gene as identified in our present study, has also been reported previously.30–32 Moreover, it was recognized as a potential cancer therapeutic target by other researchers.33 Thus, further studies are required to elucidate the functions of SGO2/CDC20 in HCC.

SKA3 (spindle and kinetochore associated complex subunit 3), also called RAMA1 or C13orf3, participates in mitosis,34 and is considered one of the genes susceptible to prostate cancer metastasis.35 Overexpression of SKA3 contributes to the progression of colorectal adenoma to carcinoma progression.36 Thus, high expression of SKA3 associated with immune infiltration in HCC suggests it is a promising candidate for prognosis evaluation and therapy of HCC.

TMCO3 (transmembrane and coiled-coil domain 3), also known as C13orf11, was identified as one of the interaction partners of EGFR (epidermal growth factor receptor).37 Considering that the EGFR-related pathway plays an important role in HCC tumorigenesis,38–40 it is necessary to clarify whether TMCO3 affects the prognosis and immune status through the EGFR-related pathway.

In a previous study, we used single-cell sequencing to analyze the landscape of infiltrating T cells in liver cancer and revealed that subsets of exhausted T cells correlated with poorer prognosis in HCC patients.41 In present study, we analyzed the relationship between the 36 key potential prognostic genes and the immune sets in HCC with TIMER. The results showed that the 36 key potential prognostic genes were significantly associated with exhausted T cell markers, including PDCD1, CTLA4, HAVCR2 and TIGIT, suggesting these genes were closely related to T cell exhaustion in TME of HCC. Previous studies demonstrated that exhausted T cells in TME highly expressed inhibitory receptors, exhibited lower effector cytokine secretion and cytolytic activity and ultimately led to failure of tumor eradication. Reversing exhausted status of T cells represents a promising breakthrough in cancer immunotherapy.42,43 Thus, the 36 genes in this study may be used as potential targets to improve the prognosis of HCC, probably through reversing T cell exhaustion in TME.

Our results also showed that the expression of these 36 genes significantly and positively correlated with PMN-MDSC (CD14-CD11b+CD33+CD15+), whose phenotypes are similar to granulocytes, whereas it negatively correlated with M-MDSC, whose phenotypes are similar to monocytes.44 Previous studies indicated that subtypes of MDSC varied in different cancer types, depending on TME.45 A lower level of MDSC in cancer was related to favorable prognosis after immune checkpoint inhibitor therapy.44,46 In a spontaneous model of uveal melanoma, PMN-MDSC infiltrated in primary tumors could promote the dissemination of cancer cells through inducing epithelial-mesenchymal transition (EMT).47 Thus, the 36 key genes may regulate the differentiation of MDSCs in TME of HCC, and influence the prognosis through up-regulating PMN-MDSC. In addition, our TIMER analysis showed that the 36 genes positively correlated with CCR8 and TGF-B1, while they negatively correlated with STAT5B. CCR8, TGF-B1 and STAT5B are well-known Treg regulating molecules, therefore these genes may affect the prognosis of HCC by participating in the differentiation of Treg in TME of HCC.

Recently, immune infiltration in TME has been extensively explored. Tumor-infiltrating lymphocytes can affect the response to immunotherapy and prognosis.48,49 A recent study indicated that immune cell infiltration within primary tumors plays a more important role than PD-L1 status in response to checkpoint inhibition in metastatic melanoma.16 Additionally, pre-existing tumor-specific T cells may have limited reinvigoration capacity through checkpoint blockade therapy;50 however, the underlying mechanisms remains unknown. In our study, the 36 key genes related to the suppressive immune status in HCC may partly elucidate their connection with the poor survival of HCC patients. To our knowledge, although prognostic value was reported in most of these 36 key genes (), their correlations with TIME in HCC were rarely explored. Better understanding of the correlation between these key prognostic genes and immune sets will open a promising avenue to ameliorate TIME and improve the immunotherapy in HCC.

Conclusions

In our study, 36 key potential prognostic genes are highly expressed in HCC and correlate with poor prognosis. These genes are significantly associated with suppressive immune status in HCC. These genes may serve as prognostic markers and potential targets to improve immunotherapy in HCC. Among them, SPINDOC, TEDC2, NUP205, SGOL2, SKA3 and TMCO3 are confirmed to be highly expressed in HCC tumors compared with adjacent normal tissues. These 6 genes may regulate proliferation of SMMC7721 (human HCC) cells. Collectively, these findings hugely enrich our knowledge of molecular characteristics and immune status in HCC.