Abstract

Introduction

Being the sixth most commonly diagnosed cancer and the second leading cause of cancer-related death, hepatocellular carcinoma (HCC) constitutes one of the major health problems worldwide [1]. There is regional heterogeneity in incidence of HCC; it is estimated that 72% of cases occur in Asia and more than 50% in China [2].

Advancements in high-throughput sequencing have significantly enhanced our comprehension of the intricate molecular landscape of HCC. Over the past decade, genome-wide association studies (GWASs) successfully identified 40 single nucleotide polymorphisms (SNPs) associated with the risk of HCC [3]. These SNPs affected several genes, including DDX18, DEPDC5, MICA, GRIK1, STAT4, KIF1B, PNPLA3, and TM6SF2 [4]. However, the functional variants and target susceptibility genes at the HCC risk regions have yet to be identified and introduced into clinical practice [4].

Up to 90% of GWAS-associated loci are distributed in intergenic regions or non-coding regions and enriched in predicted transcriptional regulatory regions, called “cis-regulatory elements (CREs)” [5]. Moreover, the GWAS-associated variants were found a high overlap with the SNPs that affect the gene expression identified in the expression quantification traits loci (eQTL) study, which suggests that many GWASs causal variants may influenced disease by altering the downstream target gene expression [5].

As one of the post-GWAS research methods, transcriptome-wide association study (TWAS) could integrate the GWAS signals and eQTLs to identified genes underlying complex traits, translating the SNP-disease associations into gene-disease associations, improving the interpretation of the genetic mechanisms underlying GWAS association signals [6]. TWAS has identified hundreds of putative susceptibility genes for many cancers, e.g. pancreatic cancer, prostate cancer, and colorectal cancer [7–10]. However, there have been relatively few TWASs published on HCC to date. [11].

Since HCC has higher incidence in East Asia than in other regions [2], we perform the TWAS in two large-scale HCC GWAS datasets from East Asian populations in this study. One is the Qidong cohort consisting of 1,161 HCC cases and 1,353 non-HCC controls who were all chronic hepatic B virus (CHB) carriers from Eastern China [12]. Another HCC GWAS dataset contains 2,122 diagnosed HCC patients and 159,201 controls from the BioBank Japan project (BBJ) [13]. To further validate our results, both single-tissue and cross-tissue TWAS have been conducted in the two East Asian populations. We hope our study will provide novel clues for the genetic mechanism of HCC in East Asia.

Patients and methods

Results

TWAS results in the Chinese population

For the Chinese Han population, our TWASs revealed 434 putative genes in the liver single-tissue analysis and 830 in the cross-tissue analysis, all with a significance level of P<0.05 (Supplementary Table 1 and Supplementary Table 2). Notably, in the liver single-tissue analysis, a single gene ACPP approached marginal significance (P_liver-tissue = 4.76×10⁻⁶; Figure 1A and Supplementary Table 1). In the cross-tissue analysis, we identified four genes significantly associated with the risk of HCC at a Bonferroni-corrected threshold of P<2.24×10⁻⁶, which were ARHGEF5 (P_cross-tissue = 2.96×10⁻⁸), AC006272.2 (P_cross-tissue = 3.61×10⁻⁷), ARNT (P_cross-tissue = 7.36×10⁻⁷), and MRPL40 (P_cross-tissue = 9.49×10⁻⁷; Figure 1B and Supplementary Table 2).

Open in a separate window

Figure 1.

Manhattan plot of the single-tissue and cross-tissue TWAS results in two East Asia populations. (A) Manhattan plot for liver single-tissue TWAS in Chinese. (B) Manhattan plot for cross-tissue TWAS in Chinese. (C) Manhattan plot for liver single-tissue TWAS in Japanese. (D) Manhattan plot for cross-tissue TWAS in Japanese. Each point represents a single gene, the x-axis represents the physical position (chromosome localization) of genes, and the y-axis represents the P-value in −log10 scale. The black solid lines represent a Bonferroni-corrected threshold of P < 2.24 × 10^-6 in (A) and (C), P < 4.14 × 10^-6 in (B) and (D); grey dotted lines represent a threshold of P < 1.00 × 10^-4. QD = Qidong GWAS cohort for Chinese Han population, BBJ = Biobank Japan GWAS cohort for Japanese, Chr = chromosome.

TWAS results in the Japanese population

For the Japanese population, TWAS identified 567 putative genes in the liver single-tissue approach and 1,342 in the cross-tissue method (P<0.05; Supplementary Table 1 and Supplementary Table 2). In the liver-tissue analysis, a single gene, ATF6B (P_liver-tissue = 9.18×10⁻⁷) was found to be significantly associated with the risk of HCC (Figure 1C and Supplementary Table 1). Additionally, the cross-tissue analysis revealed 16 genes significantly associated with the risk of HCC (Figure 1D and Supplementary Table 2).

In total, we identified 22 genes (Table 1) for which genetically predicted gene expression exhibited significant associations with HCC risk, with 21 genes passing the Bonferroni correction threshold and one coming close to it. Importantly, At three known risk loci (6p21.32, 6p21.33, and 1q42.13) identified in certain GWASs [12, 21–25], our TWAS uncovered nine genes, including ATF6B, C4B, HLA-DPA1, BTNL2, NOTCH4, HLA-DPB1, FKBPL, HLA-DRB5, and OBSCN. Whereas, the remaining 13 genes were previously unidentified.

Table 1.

Statistically significant genes selected by TWAS for HCC

					Single-tissue			Cross-tissue			Whether in new loci
GWAS cohort	Gene symbol	Region	Lead GWAS SNPa	GWAS Pa	TWAS P	FDR	TWAS Z	TWAS P	FDR	TWAS Z
Qidong	ACPP	3q22.1	rs6712250	1.73 ×10-4	4.76 ×10-6	0.057	−4.57	1.49 ×10-3	0.563	−4.57 to 4.57	Yes
	MRPL40	22q11.21	rs9606045	3.18 ×10-6	2.48 ×10-3	0.999	−3.03	9.49 ×10-7	5.29 ×10-3	−3.55 to 1.76	Yes
	ARHGEF5	7q35	rs2699503	1.48 ×10-5	3.08 ×10-5	0.185	−4.17	2.96 ×10-8	6.60 ×10-4	−4.17 to 3.54	Yes
	ARNT	1q21.3	rs75387976	2.66 ×10-3	NA	NA	NA	7.36 ×10-7	5.47 ×10-3	−2.19 to 2.19	Yes
	AC006272.2	19q13.41	rs2659048	2.04 ×10-3	NA	NA	NA	3.61 ×10-7	4.03 ×10-3	−0.54 to 2.15	Yes
BioBank Japan	ATF6B	6p21.32	rs9271377	5.29 ×10-9	9.18 ×10-7	0.011	−4.91	2.52 ×10-6	3.31 ×10-3	−5.50 to −2.69	No
	C4B	6p21.33	rs9271377	5.29 ×10-9	0.147	1.000	1.45	5.45 ×10-18	1.22 ×10-13	−3.52 to 3.83	No
	IFFO1	12p13.31	rs11829631	2.98 ×10-3	0.150	1.000	−1.44	4.90 ×10-8	1.09 ×10-4	−1.91 to 1.78	Yes
	FNDC10	1p36.33	rs181784865	5.71 ×10-4	0.304	1.000	1.03	9.54 ×10-	1.52 ×10-3	−2.25 to 3.06	Yes
	FKBPL	6p21.32	rs9271377	5.29 ×10-9	0.427	1.000	−0.79	2.45 ×10-7	4.55 ×10-4	−3.41 to 4.15	No
	LTB4R	14q12	rs28418362	8.55 ×10-6	0.571	1.000	0.57	5.62 ×10-9	2.51 ×10-5	−2.90 to 0.57	Yes
	HLA-DRB5	6p21.32	rs9271377	5.29 ×10-3	3.47 ×10-5	0.146	4.14	1.49 ×10-6	2.22 ×10-3	−2.33 to 4.56	No
	IFT57	3q13.12	rs4855552	6.09 ×10-5	4.87 ×10-4	0.346	3.49	9.16 ×10-9	3.41 ×10-5	−3.49 to 4.81	Yes
	NOTCH4	6p21.32	rs9271377	5.29 ×10-9	5.43 ×10-3	0.689	−2.78	2.47 ×10-9	1.38 ×10-5	−3.90 to 4.64	No
	CTBP1	4p16.3	rs73069980	4.11 ×10-6	6.33 ×10-3	0.697	2.73	3.34 ×10-8	1.06 ×10-4	−3.41 to 4.12	Yes
	OBSCN	1q42.13	rs115344926	8.06 ×10-7	NA	NA	NA	1.75 ×10-7	3.55 ×10-4	−3.53 to 3.35	No
	EGF	4q25	rs114251637	2.46 ×10-4	NA	NA	NA	8.42 ×10-	1.44 ×10-3	−1.11 to 2.09	Yes
	HLA-DPA1	6p21.32	rs9271377	5.29 ×10-9	NA	NA	NA	5.57 ×10-12	6.21 ×10-8	−6.02 to 2.99	No
	BTNL2	6p21.32	rs9271377	5.29 ×10-9	NA	NA	NA	1.04 ×10-9	7.73 ×10-6	−4.03 to 2.75	No
	HLA-DPB1	6p21.32	rs9271377	5.29 ×10-9	NA	NA	NA	4.10 ×10-8	1.14 ×10-4	−5.88 to 2.85	No
	LINC01405	12q24.11	rs11066015	3.80 ×10-7	NA	NA	NA	4.77 ×10-8	1.18 ×10-4	−3.68 to −2.23	Yes
	LINC01671	21q22.3	rs771148881	1.12 ×10-3	NA	NA	NA	1.60 ×10-6	2.23 ×10-3	−2.28 to 2.28	Yes

Open in a separate window

TWAS = transcriptome-wide association study, GWAS = Genome-wide association study, SNP = single nucleotide polymorphism.

All the genes listed in the above table (except ACPP) have passed Bonferroni correction threshold with a significance level of P<2.24×10⁻⁶ in cross-tissue analysis, or P<4.14×10⁻⁶ in single-tissue analysis.

^aThe lead GWAS variant and GWAS P value indicate the most statistically significant GWAS variant within ±1Mb for each gene listed.

^bGenes in reported loci implicated in previous GWASsThe black bold font: the P-value and Z-score of the genes labeled on the Manhattan plot.

Overlapped genes between Chinese and Japanese populations

Out of the 434 genes exhibiting P<0.05 in the liver single-tissue analysis for the Chinese population, 14 genes were also identified in the Japanese population with a P<0.05 in single-tissue analysis (Figure 2 and Supplementary Table 3). Furthermore, among the 830 genes with P<0.05 in the cross-tissue analysis for the Chinese population, 72 genes were likewise identified in the Japanese population with a P<0.05 in the cross-tissue analysis (Figure 2 and Supplementary Table 3). Notably, ATP6V1G2 (6p21.33), HLA-DQA1 (6p21.32), and NUAK2 (1q32.1) maintained P-values < 0.05 in both TWAS approaches and across both populations (Figure 2 and Supplementary Table 3).

Open in a separate window

Figure 2.

Venn diagram of the TWAS identified genes in different populations and approaches. Different colors represent the number of genes with P < 0.05 in different clusters. QD = Qidong cohort of Chinese, BBJ = Biobank Japan, Liver = liver single-tissue TWAS analysis (SPrediXcan), Cross = cross-tissue TWAS analysis (SMultiXcan).

Functional annotations of TWAS identified genes

In total, we identified 105 genes, comprising 22 genes with significant associations and 83 genes consistently exhibiting P<0.05 in two distinct populations, through the utilization of various TWAS methods. Of these 105 genes, 85 were categorized as protein-coding genes, which will undergo further annotation as part of our subsequent exploration of functional and molecular mechanisms.

For the functional annotation of the above 85 genes, we first carried out a pathway and process enrichment analysis with the following ontology categories using Metascape: KEGG Pathway, GO Biological Processes, GO Molecular Functions, Reactome Gene Sets, Canonical Pathways, Resource for Mammalian Protein Complex (CORUM), and Oncogenic Signatures. Genes in the top 18 significant clusters primarily engaged in immune response and transcription regulation processes, such as Staphylococcus aureus infection, initial triggering of complement, and over-expressing an important regulator of the cell cycle, Cyclin D1 (Figure 3).

Open in a separate window

Figure 3.

Visualization of the top statistically enriched terms of 85 HCC-associated genes identified by TWAS. The functional group terms were clustered based on their kappa score (≥0.3). The size of the bar represents the significance of term enrichment, and the darker the color, the more significant it is.

Subsequently, we utilized TCGA data from UALCAN to compared the mRNA expression level of 85 selected protein-coding genes in HCC tumor tissue (n=371) and normal tissues (n=50). Our analysis revealed that 50 genes, including NUAK2, ATP6V1G2, HLA-DRB5, ATF6B, FAM96B, and HSPA5, exhibited significantly higher expression in HCC tumor tissues than in normal tissues, whereas 11 genes, such as CAT, PMAIP1, and RXRG, displayed significantly lower expression in HCC tumor tissue than in normal tissues (P<0.01; Table 2).

Table 2.

Functional examination of the TWAS-identified genes associated with HCC

Gene symbol	Region	Whether in new loci	Meta-analysis		Differential expression analysis		Survival analysis			Number of gene effect score < −0.5
			TWAS Z	TWAS P	P value	Expression in HCC compared to controla	OS; HR (P)	DFS; HR (P)	High expression and survival timeb	CRISPR	RNAi
NUAK2	1q32.1	Yes	2.86	4.23×10−3	4.50×10−6	Up	1.50 (0.035)	1.20 (0.160)	Shorter	0	0
HLA-DQA1	6p21.32	No	3.27	1.07×10−3	7.82×10−2	Down	0.75 (0.110)	0.74 (0.053)	NA	0	/
ATP6V1G2	6p21.33	No	−3.27	1.06×10−3	4.67×10−10	Up	1.10 (0.700)	0.95 (0.760)	NA	0	1
HLA-DRB5	6p21.32	No	3.01	2.58×10−3	2.30×10−3	Up	0.84 (0.320)	0.84 (0.240)	NA	0	0
HLA-DRA	6p21.32	No	1.65	9.96×10−2	2.09×10−2	Down	0.93 (0.690)	0076 (0.079)	NA	0	0
FAM96B	16q22.1	Yes	3.18	1.47×10−3	1.62×10−12	Up	1.30 (0.120)	1.40 (0.018)	Shorter	20	0
UGT2B17	4q13.2	Yes	−1.61	0.110	NA	NA	0.90 (0.560)	1.10 (0.660)	NA	0	0
VTI1A	10q25.2	No	−2.58	9.75×10−3	<1×10−12	Up	1.20 (0.420)	1.20 (0.280)	NA	0	0
LTBP3	11q13.1	Yes	−1.62	0.110	1.43×10−9	Up	1.10 (0.680)	0.80 (0.160)	NA	1	0
SNX27	1q21.3	Yes	−2.21	0.027	<1×10−12	Up	1.30 (0.140)	1.20 (0.310)	NA	0	0
KLHL12	1q32.1	Yes	−2.07	0.038	<1×10−12	Up	1.40 (0.055)	1.10 (0.540)	NA	0	0
HSPA5	9q33.3	Yes	0.43	0.660	<1×10−12	Up	1.50 (0.017)	1.10 (0.680)	Shorter	20	1
AGPAT1	6p21.32	No	−1.88	0.061	1.62×10−12	Up	1.30 (0.130)	1.30 (0.120)	NA	0	0
PRMT7	16q22.1	Yes	−0.59	0.560	1.62×10−12	Up	1.00 (0.810)	1.10 (0.370)	NA	1	1
PKD1L3	16q22.2	Yes	2.28	0.023	8.31×10−1	Down	0.61 (0.020)	0.66 (0.025)	Longer	/	0
CAT	11p13	Yes	2.79	5.32×10−3	1.90×10−12	Down	0.63 (0.011)	0.84 (0.250)	Longer	0	0
DUSP3	17q21.31	Yes	−2.11	0.035	<1×10−12	Up	1.20 (0.250)	1.30 (0.120)	NA	0	0
SOCS6	18q22.2	Yes	−2.50	0.013	4.12×10−3	Down	1.40 (0.051)	1.20 (0.280)	NA	0	0
ACPP	3q22.1	Yes	−1.73	0.084	2.61×10−4	Down	0.97 (0.870)	0.87 (0.400)	NA	0	0
MRPL40	22q11.21	Yes	−1.47	0.141	1.62×10−12	Up	1.10 (0.690)	1.30 (0.140)	NA	5	4
ARHGEF5	7q35	Yes	−2.30	0.021	1.90×10−11	Up	1.20 (0.430)	1.10 (0.400)	NA	4	1
ARNT	1q21.3	Yes	−1.80	0.072	1.62×10−12	Up	1.30 (0.150)	1.10 (0.580)	NA	0	0
ATF6B	6p21.32	No	−3.94	8.19×10−5	<1×10−12	Up	1.40 (0.060)	1.20 (0.230)	NA	0	0
C4B	6p21.33	No	4.81	1.49×10−6	NA	NA	0.85 (0.360)	0.84 (0.270)	NA	/	0
IFFO1	12p13.31	Yes	2.77	5.55×10−3	<1×10−12	Up	0.89 (0.520)	1.00 (0.750)	NA	0	0
FNDC10	1p36.33	Yes	2.91	3.56×10−3	NA	NA	NA	NA	NA	1	/
FKBPL	6p21.32	No	−3.10	1.93×10−3	<1×10−12	Up	1.50 (0.018)	1.70 (0.001)	Shorter	0	0
LTB4R	14q12	Yes	−3.53	4.23×10−4	1.62×10−12	Up	1.10 (0.450)	1.30 (0.056)	NA	2	0
HLA-DRB5	6p21.32	No	3.01	2.58×10−3	2.30×10−3	Up	0.84 (0.320)	0.84 (0.240)	NA	/	/
IFT57	3q13.12	Yes	3.40	0.001	7.13×10−8	Up	1.10 (0.470)	0.97 (0.850)	NA	1	0
NOTCH4	6p21.32	No	3.43	6.06×10−4	1.11×10−16	Up	0.90 (0.550)	1.20 (0.170)	NA	0	0
CTBP1	4p16.3	Yes	3.53	4.12×10−4	<1×10−12	Up	1.20 (0.390)	1.00 (0.890)	NA	0	1
OBSCN	1q42.13	No	−3.07	2.13×10−3	1.86×10−12	Up	1.50 (0.015)	1.10 (0.680)	Shorter	0	0
EGF	4q25	Yes	−3.30	9.54×10−4	4.72×10−12	Up	1.60 (0.015)	0.93 (0.630)	Shorter	0	0
HLA-DPA1	6p21.32	No	−4.06	4.82×10−5	1.97×10−1	Down	0.93 (0.670)	0.80 (0.150)	NA	0	0
BTNL2	6p21.32	No	−3.95	7.96×10−5	NA	NA	NA	NA	NA	/	0
HLA-DPB1	6p21.32	No	−2.89	3.84×10−3	8.72×10−2	Down	1.00 (0.930)	0.80 (0.150)	NA	0	0
PHYH	10p13	Yes	−2.07	0.038	4.39×10−4	Down	0.96 (0.820)	1.10 (0.730)	NA	0	0
HABP2	10q25.3	Yes	−0.80	0.423	2.85×10−6	Down	0.86 (0.410)	0.89 (0.440)	NA	0	0
SSSCA1	11q13.1	Yes	−2.58	9.81×10−3	1.62×10−12	Up	1.10 (0.520)	1.10 (0.630)	NA	/	/
ARHGEF17	11q13.4	Yes	−1.02	0.306	1.79×10−12	Up	1.20 (0.240)	0.92 (0.580)	NA	0	/
H3F3C	12p11.21	Yes	−2.17	0.030	1.40×10−4	Up	0.90 (0.570)	0.96 (0.820)	NA	/	/
CSRNP2	12q13.12	Yes	2.30	0.022	1.62×10−12	Up	1.60 (0.008)	1.30 (0.065)	Shorter	0	/
TFCP2	12q13.12-q13.13	Yes	2.41	0.016	1.62×10−12	Up	1.40 (0.074)	1.20 (0.330)	NA	0	/
CCDC38	12q23.1	Yes	−1.13	0.260	1.16×10−1	Down	0.68 (0.030)	0.81 (0.160)	Longer	0	0
UNG	12q24.11	Yes	1.98	0.047	<1×10−12	Up	1.40 (0.082)	1.70 (0.001)	Shorter	0	0
GP1BA	17p13.2	Yes	2.53	0.011	1.99×10−1	Down	0.73 (0.080)	0.70 (0.023)	Longer	0	0
KRT24	17q21.2	Yes	−0.70	0.483	NA	NA	NA	NA	NA	0	1
SLC25A39	17q21.31	Yes	−1.04	0.298	1.62×10−12	Up	1.90 (0.001)	1.50 (0.015)	Shorter	0	0
EPX	17q22	Yes	1.34	0.179	5.66×10−2	NA	NA	NA	NA	0	0
STXBP4	17q22	Yes	−0.61	0.543	1.62×10−12	Up	1.50 (0.018)	1.40 (0.018)	Shorter	0	0
PMAIP1	18q21.32	Yes	1.96	0.050	1.32×10−3	Down	1.50 (0.020)	1.30 (0.120)	Shorter	0	0
XAB2	19p13.2	Yes	−1.00	0.315	1.62×10−12	Up	1.10 (0.740)	1.20 (0.350)	NA	20	16
HPCAL4	1p34.2	Yes	−2.12	0.034	2.19×10−1	Down	NA	NA	NA	0	0
C1QC	1p36.12	Yes	−2.14	0.033	7.74×10−2	Down	1.20 (0.320)	0.86 (0.340)	NA	0	0
WFDC13	20q13.11	Yes	−2.12	0.034	3.20×10−1	NA	NA	NA	NA	0	0
SHROOM3	4q21.1	Yes	−3.26	1.10×10−3	3.12×10−3	Down	0.95 (0.770)	1.10 (0.720)	NA	0	0
NR3C1	5q31.3	Yes	−2.39	0.017	9.78×10−5	Up	1.10 (0.780)	1.30 (0.090)	NA	0	0
PCDHB1	5q31.3	Yes	0.57	0.570	3.70×10−2	NA	NA	NA	NA	0	0
HCG20	6p21	Yes	−0.90	0.369	NA	NA	1.50 (0.021)	1.10 (0.480)	Shorter	/	/
PBX2	6p21.32	No	−3.10	1.94×10−3	1.62×10−12	Up	1.10 (0.460)	1.10 (0.480)	NA	1	0
EGFL8	6p21.32	No	1.95	0.051	<1×10−12	Up	0.85 (0.370)	1.40 (0.031)	Shorter	0	0
HLA-DRB1	6p21.32	No	1.65	0.099	6.62×10−4	Up	0.84 (0.310)	0.87 (0.380)	NA	0	0
CFB	6p21.33	No	1.93	0.054	4.59×10−2	Down	0.66 (0.020)	0.79 (0.130)	Longer	0	0
CYP21A2	6p21.33	No	−3.25	1.17×10−3	5.33×10−13	Up	1.50 (0.037)	1.10 (0.670)	Shorter	0	1
DDR1	6p21.33	No	2.35	0.019	1.02×10−9	Up	1.20 (0.220)	1.10 (0.530)	NA	0	0
MICB	6p21.33	No	−2.59	9.51×10−3	1.62×10−12	Up	1.30 (0.130)	1.60 (0.002)	Shorter	0	0
MRPS18B	6p21.33	No	1.29	0.199	<1×10−12	Up	1.20 (0.250)	1.00 (1.000)	NA	10	0
NFKBIL1	6p21.33	No	1.08	0.281	<1×10−12	Up	1.20 (0.270)	1.40 (0.470)	NA	0	0
SIM1	6q16.3	Yes	2.32	0.021	2.66×10−3	Down	0.91 (0.600)	0.91 (0.570)	NA	0	0
GRB10	7p12.1	Yes	1.35	0.175	<1×10−12	Up	1.10 (0.610)	0.94 (0.700)	NA	0	0
GRID2IP	7p22.1	Yes	−0.23	0.822	5.83×10−4	Down	0.90 (0.640)	0.93 (0.700)	NA	0	/
DEFA4	8p23.1	Yes	−2.06	0.040	7.69×10−1	NA	NA	NA	NA	0	0
FAM86B2	8p23.1	Yes	−2.38	0.017	4.24×10−3	Down	1.10 (0.590)	1.30 (0.085)	NA	0	/
RABEPK	9q33.3	Yes	0.38	0.707	2.74×10−14	Up	0.91 (0.600)	1.00 (0.830)	NA	0	0
CCDC62	12q24.31	Yes	3.32	8.97×10−4	4.62×10−4	Up	1.20 (0.390)	1.20 (0.260)	NA	0	0
HCAR2	12q24.31	Yes	3.16	1.57×10−3	7.90×10−2	Down	1.20 (0.350)	0.94 (0.710)	NA	0	0
ITGAD	16p11.2	Yes	1.14	0.254	6.20×10−1	Down	0.70 (0.045)	0.75 (0.061)	Longer	0	0
PRRT2	16p11.2	Yes	−0.74	0.461	<1×10−12	Up	1.10 (0.530)	1.20 (0.240)	NA	0	0
POLRMT	19p13.3	Yes	2.73	6.33×10−3	1.62×10−12	Up	1.30 (0.120)	1.40 (0.032)	Shorter	9	/
RXRG	1q23.3	Yes	1.04	0.300	2.31×10−5	Down	0.77 (0.150)	0.68 (0.013)	Longer	0	1
MPHOSPH10	2p13.3	Yes	0.71	0.480	<1×10−12	Up	1.20 (0.300)	1.40 (0.027)	Shorter	20	6
SH3YL1	2p25.3	Yes	2.74	6.24×10−3	8.80×10−2	Down	1.00 (0.830)	0.87 (0.380)	NA	0	0
RABL2A	2q14.1	Yes	3.05	2.27×10−3	1.62×10−12	Up	1.40 (0.079)	1.60 (0.004)	Shorter	4	1
GRPEL1	4p16.1	Yes	1.02	0.306	1.90×10−1	Down	1.00 (1.000)	0.83 (0.220)	NA	20	6

Open in a separate window

TWAS = transcriptome-wide association study, OS = overall survival, DFS = disease-free survival, HR (P) = hazard ratio (P-value), RNAi = RNA interference.

^a“Up” means that the expression levels of these genes were significantly up-regulated in hepatocellular carcinoma tissues compared to normal tissues; “Down” means the genes were down-regulated in HCC tissues.

^b“shorter” means that the patients with high mRNA expression levels of these genes had shorter survival time compared with those with low expression level; “longer” means that patients with high expression levels of genes had longer survival time.

We then employed GEPIA2 tool to perform survival analyses for the above 85 genes. In the overall survival (OS) analyses, we found that patients with high expression levels of NUAK2, HSPA5, FKBPL, OBSCN, EGF, CSRNP2, SLC25A39, STXBP4, PMAIP1, HCG20, and CYP21A2 had worse prognoses than those with low expression levels of these 11 genes. Conversely, patients with high expression levels of PKD1L3, CAT, CCDC38, CFB, and ITGAD had better prognoses than those with low expression levels of these five genes (Table 2; P<0.05). In the disease-free survival (DFS) analyses, patients with low expression of 10 genes (FAM96B, FKBPL, UNG, SLC25A39, STXBP4, EGFL8, MICB, POLRMT, MPHOSPH10, and RABL2A) and high expression of three genes (PKD1L3, GP1BA, and RXRG) exhibited better prognoses compared with other patients (Table 2).

Finally, we assessed the dependency of the 85 genes in 23 HCC cell lines for both CRISPR-Cas9 and RNAi data in the DepMap database (Supplementary Table 4). We employed the gene effect scores to determine the gene’s dependency on a specific cell line or tissue. Negative gene effect scores indicated that the cell line’s growth would slow down after experimental manipulation of a specific gene, with more negative scores indicating stronger dependency on the cell line. We set the cut-off value at −0.5 according to the DepMap website (https://forum.depmap.org/t/depmap-genetic-dependencies-faq/131). The results revealed that for FAM96B, HSPA5, XAB2, MPHOSPH10, and GRPEL, 20 out of 23 (86.96%) cell lines fell below the cut-off value (−0.5) for the gene effect score in CRISPR-Cas9 data (Table 2). This suggests that HCC cells exhibited a high sensitivity to these five genes, and the growth of cancer cells was significantly constrained following CRISPR-Cas9 editing, highlighting their potential as promising drug targets for HCC.

Discussion

Utilizing large-scale GWAS data of our cohort of Chinese CHB patients with HCC and BBJ HCC consortium, we conducted the TWAS analysis in East Asian population to explore novel susceptibility loci and putative causally relevant genes for HCC risk. In our current study, we identified 22 putative susceptibility genes for HCC after multiple testing corrections. Notably, 13 of these genes were located in the regions previously unexplored in GWASs. Furthermore, we conducted an overlapped analysis of TWAS-identified genes with a significance level of P<0.05. Ultimately, 83 genes were identified in both Chinese and Japanese populations. Moreover, our functional annotations revealed some potentially pathogenic genes that may serve as underlying risk factors for HCC. Specifically, FAM96B at 16q22.1, HSPA5 at 9q33.3, POLRMT at 19p13.3, MPHOSPH10 at 2p13.3, and RABL2A at 2q14.1 displayed consistent effect directions in various functional aspects.

Our findings provide further support for several genes that have been previously implicated in GWASs or other researches. Earlier studies have identified 31 susceptibility loci associated with hepatitis virus-related HCC, nine of which were located in the major histocompatibility complex (MHC) region [12, 21–24, 26], suggesting the pivotal role of the anti-infectious immunity in the progression of HCC. In our TWAS results, we reported that eight of the 22 significant genes (ATF6B, C4B, HLA-DPA1, BTNL2, NOTCH4, HLA-DPB1, FKBPL, and HLA-DRB5) were situated within the MHC region, indicating a potential impact on the pathogenesis of HCC through the immune processes related to viral infection. Additionally, six of the 22 genes (ARNT, LTB4R, IFT57, CTBP1, OBSCN, EGF) have previously been shown to be associated with HCC through bioinformatics methods or experimental studies [27–33].

Through a comprehensive TWAS across diverse populations and using different TWAS analysis methods, we identified 83 susceptible genes that exhibited overlap in at least one TWAS approach within the two populations. Among these, compelling evidence supported the association of three genes (NUAK2, HLA-DQA1, and ATP6V1G2) with the risk of HCC, as they were consistently validated in both two populations using distinct TWAS analysis approaches. Notably, variants such as rs9272105 in HLA-DQA1 have been proven to be associated with the susceptibility of HCC in previous report of GWAS [21]. NUAK2 functions as a stress response kinase involved in multiple processes, including cellular response to glucose starvation, negative regulation of apoptotic process, and protein phosphorylation [34]. In HCC, it serves as an essential mediator of liver proliferation and tumorigenesis by targeting Hippo/YAP signaling pathway [35]. Known as ATPase H+ Transporting V1 Subunit G2, ATP6V1G2 encodes an enzyme transporter V-ATPase, which regulates the acidification of intracellular compartments in eukaryotic cells [36]. A previous study indicated that high expression of the ATPase-V0a3 subunit (TCIRG1) was observed in HCC and promoted tumor cell growth and proliferation [37, 38]. Although ATP6V1G2 has been proven to be a risk gene for colorectal cancer [39] and glioma [38], there is currently no evidence linking ATP6V1G2 to HCC progression.

In our analysis, we found C4B (6p21.33) may be the potential candidate risk gene in HCC. C4B encodes the basic form of complement factor 4, and together with the C4A gene, is part of the classical activation pathway [40]. C4A/C4B up-regulation is thought to contribute to HCC development through inflammatory and immunosuppressive mechanisms [41]. In the context of HCC, the specific role of the C4B gene may be related to its involvement in the immune response against cancer cells. But the exact role of C4B in HCC can vary and may be influenced by various factors, including the specific characteristics of the tumor and the individual’s immune response. We need further research to fully understand the precise mechanisms by which genes may influence the HCC.

The pathway enrichment analysis revealed several candidate pathways for HCC, some of which have been reported to be directly or indirectly involved in HCC development. Notably, the cell cycle, immune response, complement reactions, metabolic processes, transcription factor activity, and the oxidative stress pathways have been widely recognized as significant contributors to HCC [42–45]. Additionally, researchers have provided evidence of the enrichment of the Staphylococcus aureus infection pathway in HCC tumor samples [46].

Functional annotation revealed that five genes (FAM96B, HSPA5, POLRMT, MPHOSPH10, and RABL2A) consistently displayed significant effects across various methods. For instance, their expression was significantly higher in HCC tissues than in normal controls. Patients with high expression of these genes experienced shorter survival times and a poorer prognosis, compared with others. Furthermore, most of the HCC cell lines demonstrated sensitivity to these five genes. And their knockout led to a restriction in the growth of HCC cells, highlighting the potential role of these five genes as oncogenes (Table 2). FAM96B, also called cytosolic iron-sulfur assembly component 2B (CIAO2B), is a highly conserved homologous protein belonging to the MIP18 family. An important paralog of FAM96B is FAM96A/CIAO2A. FAM96A and FAM96B participate in regulating tumor cell apoptosis by interacting with different proteins [47]. Some studies have shown that FAM96A can act as a tumor suppressor gene in various human cancers, including liver cancer [48]. But there is currently no detailed research report on the functional role of FAM96B in tumor progression at present. In one Chinese research, FAM96B expression was found to be down-regulated in HCC [47]. But in our study, we observed an up-regulation of FAM96B expression in TCGA-LIHC datasets. Further In vitro/in vivo experiments are needed to elucidate the molecular regulatory mechanisms and functional roles of FAM96B in the occurrence and development of HCC in different populations. Located in the lumen of the endoplasmic reticulum, HSPA5 is a member of the heat shock protein 70 (HSP70) family, and it is involved in endoplasmic reticulum protein folding and assembly. HSPA5 plays a role in endoplasmic reticulum stress and is implicated in various human diseases, including virus infection and tumorigenesis [49, 50]. In HCC, HSPA5 interacts with CD5L, a soluble scavenger cysteine-rich protein that modulates inflammatory responses, and promotes the proliferation of liver cancer cells, playing an antiapoptotic role [51]. In both Hispanic and Asian populations, HSPA5 could be the potential diagnostic biomarkers for HCC patients [52]. As an essential primer synthesis to initiate mitochondrial DNA (mtDNA) replication, POLRMT play a key role in the replication of mammalian mtDNA and is part of a transcriptional mechanism that controls the switch between primer formation for mtDNA replication and mitochondrial gene expression [53]. However, there is currently no direct evidence or functional experiments conducted on HCC samples. The role of POLRMT in HCC deserves further investigation. MPHOSPH10 is a M-Phase Phosphoprotein involved in rRNA processing in the nucleus and cytosol. Several GWASs for neuropsychiatric diseases have identified several SNPs on MPHOSPH10 [54, 55]. But the role of MPHOSPH10 in HCC has not been experimentally confirmed and requires further exploration. The RABL2A gene is a member of the RAS oncogene family-like 2 and plays an essential role in the proper function of flagella and cilia [56]. For HCC, one study illustrated that the expression of RABL2A was significantly up-regulated in sorafenib-resistant HCC cells, suggesting the involvement of the RABL2A–CCDC34 axis in mediating p38/MAPK and JNK/MAPK signaling pathways, thereby contributing to acquired sorafenib resistance in HCC [57].

It should be noted that there are several limitations in this study. Firstly, of the 22 genes with significant threshold values, the genes identified in the Chinese population were not validated in the Japanese cohort, even with a nominal significance. This discrepancy may be attributed the differing primary causes of HCC in the two populations. The main cause of HCC in China is hepatitis B viral (HBV) infection, while in Japan it is the hepatitis C viral (HCV) infection [58, 59]. Secondly, our HCC GWAS data were both derived for Asian populations, but the TWAS reference weights were based on European populations. To validate our results, it is essential to establish TWAS reference weights specifically tailored to Asian populations. Thirdly, though the TWAS results in Chinese and Japanese found three overlapping genes using different methods, these associations were only nominally significant and may include false positives. Further investigation into the functions of these genes in HCC is warranted.

In summary, we conducted a comprehensive TWAS in diverse populations using different TWAS approaches, identified some novel candidate genes and pathways associated with HCC. These findings contribute to a better understanding of the genetic mechanisms underlying HCC. But still need further biological studies to confirm our findings in the future.

Conclusion

TWASs have the potential ability to uncover new candidate genes that go beyond the associations identified in GWASs. Our research offers a valuable genetic insight into the biology of HCC.

Supplementary Material

Acknowledgements

Contributor Information

Jingjing Zhang, School of Public Health (Shenzhen), Sun Yat-sen University, Guangzhou, Guangdong, P. R. China. School of Public Health (Shenzhen), Shenzhen campus of Sun Yat-sen University, Shenzhen, Guangdong, P. R. China.

Qingrong Zhang, School of Public Health (Shenzhen), Sun Yat-sen University, Guangzhou, Guangdong, P. R. China. School of Public Health (Shenzhen), Shenzhen campus of Sun Yat-sen University, Shenzhen, Guangdong, P. R. China.

Wenyan Hu, School of Public Health (Shenzhen), Sun Yat-sen University, Guangzhou, Guangdong, P. R. China. School of Public Health (Shenzhen), Shenzhen campus of Sun Yat-sen University, Shenzhen, Guangdong, P. R. China.

Yuxuan Liang, School of Public Health (Shenzhen), Sun Yat-sen University, Guangzhou, Guangdong, P. R. China. School of Public Health (Shenzhen), Shenzhen campus of Sun Yat-sen University, Shenzhen, Guangdong, P. R. China.

Deke Jiang, Department of Infectious Diseases and Hepatology Unit, Nanfang Hospital, Southern Medical University, Guangzhou, Guangdong, P. R. China.

Haitao Chen, School of Public Health (Shenzhen), Sun Yat-sen University, Guangzhou, Guangdong, P. R. China. School of Public Health (Shenzhen), Shenzhen campus of Sun Yat-sen University, Shenzhen, Guangdong, P. R. China.

References