The transcriptomic data and corresponding clinical information of HCC cases were acquired from TCGA liver hepatocyte carcinoma (TCGA-LIHC) project till April 26, 2021, and were used as the training set. It included 365 tumorous specimens and 50 adjacent non-tumor tissues ( Supplementary Table 1 ). The transcriptomic data and clinical information of 240 HCC cases were obtained from the International Cancer Genome Consortium (ICGC) portal (LIRI-JP project, https://dcc.icgc.org), and was used as a validation set ( Supplementary Table 1 ). The “limma” package in R was used to analyze the mRNA expression profile of MED8 in tumor and adjacent non-tumor tissues in TCGA-LIHC cohort. The expression of MED8 among 50 paired tumor tissue and the corresponding adjacent non-tumor specimens was compared. The correlation of MED8 expression with the clinicopathological features was explored. TCGA-LIHC cohort cases were stratified into low- or high- groups based on the median expression of MED8, the difference in the overall survival (OS) between the groups was analyzed. Univariate and multivariate Cox proportional hazards regression analyses were used to identify the independent factors that affecting OS in TCGA-LIHC patients.

TCGA-LIHC cohort was stratified into high- or low- groups based on the median mRNA expression of MED8. The potential pathways upregulated due to high MED8 expression were identified by GSEA. The annotated gene sets from c2.cp.kegg.v7.4.symbols.gmt were used as the reference. The number of gene set permutations was selected as 1,000 in each analysis. The significant pathways were filtered by the normalized enrichment scores (NES) > 1, false discovery rate < 0.25 and nominal p < 0.05. The CIBERSORT analysis is a deconvolution algorithm that predicts the 22 immune cell type proportions from gene expression signature (16). The proportions of tumor-infiltrating immune cells in high- and low- groups was inferred by the CIBERSORT online tool (https://cibersort.stanford.edu/runcibersort.php), followed by filtration with the criterion of P < 0.05. The correlation between the level of MED8 expression and the proportion of tumor-infiltrating immune cells was examined using the Spearman correlation test.

The candidate MED8-related immunomodulators including immunoinhibitors and immunostimulators were retrieved from TISIDB (http://cis.hku.hk/TISIDB/index.php), a website for predicting immune-tumor interactions (17). The inclusion criteria for the significantly related immunomodulators were based on the Spearman correlation results and P-value < 0.05. First, the univariate analysis was used along with the Cox proportional hazards regression model to identify significant immunomodulators that affected the OS of patients in TCGA-LIHC cohort. The resulting immunomodulators filtered using univariate analysis were used as the input for multivariate Cox analysis. Thus, a predictive model was constructed and the stepwise variables were selected with the Akaike information criteria executed using the “stepAIC” algorithm in the “MASS” package of R (18). The final prediction model was built based on the optimal variables. Risk scores considered as the prognostic index and were calculated from the equation: risk score = gene 1 expression × gene 1 coefficient + gene 2 expression × gene 2 coefficient + … + gene n expression × gene n coefficient. The coefficient of each gene was obtained from the Cox model. The cut-off for the risk scores came from the predictive model that separated the cases into low-risk and high-risk groups and was identified in X-tile software (version 3.6.1).

The clustering abilities of the risk scores derived from TCGA-LIHC cohort were visualized using the t-distributed stochastic neighbor embedding (t-SNE) algorithm and principal component analysis. Survival differences between low-risk and high-risk groups were assessed using Kaplan-Meier curves and log-rank test. Stratified OS analyses were conducted to assess the differences in the subgroups of early clinicopathology classifications including T1 stage, N0 stage, M0 stage, pathological stage 1, and histological grade 1. The prediction accuracy of the predictive model in 1-year, 3-year, and 5-year OS was assessed and compared with the TNM model that incorporated stage and grade by measuring the area under the curve (AUC) of the receiver operating characteristic (ROC) curve via the “timeROC” package of R. ROC curve analyses were used to determine whether our model had better prediction accuracy relative to the previously published immune-related models for prediction of long-term survival. The “Meng model”, “Dai model”, and “Wang model” were included in the comparison, which were reported to have a good performance in prognostic prediction (7, 10, 11).

A validation set from the ICGC cohort was adopted to externally validate the prognostic performance of the model. The risk scores were figured for each patient in the ICGC cohort with the same formula, and the patients were then categorized into low-risk and high-risk groups using the same cut-off. KM curve, ROC, PCA, and t-SNE analyses was performed as described above. ROC analyses were only conducted for 1-year and 3-year OS due to the small number of cases with the OS time of more than 5-year in the ICGC cohort.

To verify the predictive value of the model for recurrence of HCC, the TCGA-LIHC cohort specimens were classified into low-risk and high-risk groups in line with the above cut-off. Significant differences for comparing the disease-free survival (DFS) time were assessed with by the KM curve analysis. The performance of the predictive model in DFS was estimated using ROC analysis as described above.

For quantifying the survival probability for each patient in TCGA-LIHC cohort, the nomogram was developed based on TCGA-LIHC cohort by incorporating the clinicopathological features and risk scores of patients using the R “rms” package. Plotting a calibration described the consistency between the predicted and actual probabilities of the OS through the bootstrap method with 1,000 replicates.

The power of the predictive model in predicting ICB therapeutic responses in HCC patients was inferred from the immunophenoscore (IPS) derived from the Cancer Immunome Atlas (https://tcia.at/home) and the tumor immune dysfunction and exclusion (TIDE) score derived from the TIDE online website (http://tide.dfci.harvard.edu/). The IPS value positively reflected the tumor immunity from the expression of crucial components of tumor immunogenicity; the TIDE prediction scores represented the potential for tumor immune escape. The lower the TIDE score and the higher the IPS value, the better the therapeutic response to ICB. Single sample gene set enrichment analysis (ssGSEA) was used to examine the diversities in the activities of 13 immune-related pathways and infiltration scores for 16 immune cells types between the high-risk and low-risk groups of TCGA-LIHC cohort using the R “GSVA” package. The predictive capacity of the model for response to chemotherapeutic and targeted drugs was evaluated using R “pRRophetic” package and the Genomics of Drug Sensitivity in Cancer web tool (https://www.cancerrxgene.org/) (19). The half-maximal inhibitory concentration (IC50) represented the chemotherapeutic sensitivity in the different risk groups of TCGA-LIHC. To examine the differences in the biological pathways between the low-risk and high-risk groups of TCGA-LIHC cohort, GSEA was conducted as described above.

Furthermore, the ICGC cohort was also adopted to externally verify the predictive capacity of the model in responding therapeutic sensitivity. The ICB therapy response of low-risk and high-risk groups of ICGC cohort was only inferred by the TIDE value due to the restricted dataset of the Cancer Immunome Atlas.

The prepared paraffin-embedded sections of HCC and non-tumor liver tissues were sequentially deparaffinized, rehydrated, and antigen retrieval was performed by by heating treatment in citrate buffer. After blocking with 3% bovine serum albumin, the sections were stained with the primary antibody, and polyclonal rabbit anti-MED8 (1:200, Affinity), and kept overnight at 4°C. The following day after the primary antibody incubation, the sections were washed and stained with horseradish peroxidase-conjugated secondary antibody (1:200, Abcam) for one hour at 37°C and visualized using diaminobenzidine.

Human HCC cell lines, HepG2, and Huh7 were obtained from the cell bank of the Chinese Academy of Science (Shanghai, China). For the generation of HepG2 with MED8 knocked down (HepG2 KD) and Huh7 with MED8 knocked down (Huh7 KD), cells were transfected with 50nM specific siRNA (Gene-Pharma, China) constructs for human MED8 sequence using Lipofectamine 2000 (Invitrogen, USA). The scrambled siRNA construct was used as the negative control (NC). Cells were cultured in Dulbecco’s modified eagle medium supplemented with 10% fetal bovine serum (37°C and 5% CO₂) and transfected when the confluency reached 50%. The specific siRNA constructs for MED8 are listed in Supplementary Table 2 .

Total RNA from cells was isolated by the TRIzol method (Invitrogen, USA) and dissolved in 50μl RNase-free water. The RNA (1 ug) was reverse-transcribed using the reverse transcription kit (Takara, Japan). The qRT-PCR reaction was performed using the SYBR Premix Ex Taq kit (Takara, Japan) according to the manufacturer’s instructions. The relative differences in mRNA expression were measured by the 2^−ΔΔCT method. Primer sequences of the genes adopted for qRT-PCR are specified in Supplementary Table 3 .

The total protein quantity from HCC cells was determined by the bicinchoninic acid assay. The proteins were separated on sodium dodecyl sulfate gel and transferred onto the PVDF membranes. Proteins on membranes were blocked with 5% non-fat milk for 2 hours, probed with primary antibodies (1: 1000) overnight at 4°C and thereafter incubated with the secondary antibody (1: 10000) for 1 hour. The enhanced chemiluminescence (Beyotime) was used to visualize the brands. The antibodies used are specified in Supplementary Table 4 .

Cells of HepG2 KD, HepG2 NC, Huh7 KD and Huh7 NC groups were harvested and seeded in 96-well plates at a density of 2,000 cells per well. The 10μl of cell counting kit-8 solution (Yeasen) was added to the cells at 12, 24, 36, 48, 72, and 96 hours after seeding before measuring the absorbance at 450nm. For the clonogenic assay, cells of each group were seeded in six-well plates at 1,000 cells per well. After incubation for two weeks, cells were fixed in 4% paraformaldehyde for 30 minutes, stained with 0.1% crystal violet solution, and counted.

Cell migration abilities of cells treated with specific siRNA or scramble siRNA constructs were assessed by the Transwell chambers (Corning) assays. A total of 1 × 10^5 cells were seeded onto the upper chamber. After incubation for 48 hours, the migrated cells were fixed in 4% paraformaldehyde and stained with crystal violet solution.

All data analyses were conducted using the R software (Version 4.0.3) and Prism software (GraphPad software Inc., version 9.0.0). The data are presented as mean ± standard deviation. In vitro assays were independently repeated a minimum of three times. Pairwise comparisons were performed using Student’s t-test for normally distributed variables or the Wilcoxon rank-sum test for non-normally distributed variables. P-values < 0.05 were considered statistically significant.

The mRNA expression of MED8 was significantly elevated in TCGA-LIHC tissues relative to the adjacent non-tumor liver tissues ( Figures 1A, B ). According to the HPA database ( Figures 1C, D ) and IHC staining ( Figures 1E, F ), the MED8 protein expression was enhanced in HCC. Additionally, the high expression of MED8 was associated with poor OS ( Figure 1J ) and the advanced histological grade ( Figures 1G –I ) of TCGA-LIHC cohort. Univariate and multivariate identified high expression of MED8 as a detrimental independent prognostic factor for OS in HCC patients ( Figures 1K, L ).

In the aforementioned results, high MED8 expression was significantly correlated with poor prognosis and advanced TNM parameters in HCC. GSEA was used to examine the pathways upregulated due to high MED8 expression in HCC. The results of GSEA showed that the top three enrichment pathways according to the NES scores were in the processes of the cell cycle, oocyte meiosis, and spliceosome ( Figure 1M ). The results of CIBERSORT analysis showed that proportions of naïve B cells and resting memory CD4+ T cells was negatively correlated with MED8 expression and those of activated memory CD4+ T cells and follicular helper T cells exhibited positive correlations in HCC ( Figure 1N ).

To characterize the pathophysiological significance of MED8 in HCC cell lines, both mRNA and protein expressions of MED8 were suppressed in HepG2 cells ( Figures 2A –C ) and Huh7 cells ( Figures 2E –G ) using siRNA constructs. Proliferating cell nuclear antigen (PCNA), a cell proliferation marker, reflected the status of cell proliferation. The protein level of PCNA reduced significantly after siRNA treatment. The results of CCK-8 assay indicated that the MED8 knockdown inhibited the proliferation rates in HepG2 cells and Huh7 cells compared to their respective controls ( Figures 2D, H ). Furthermore, siRNA-mediated MED8 knockdown could significantly inhibit the clonogenicity ( Figures 2I, J ) and migration ( Figures 2K, L ) in HepG2 cells and Huh7 cells. These results indicated that MED8, indeed, played an essential role in the maintenance of proliferation and migration in the HCC cells, which suggested its potential as a promising target for HCC treatment.

In total, 18 immunostimulators ( Supplementary Figure S1 ) and 7 immunoinhibitors ( Supplementary Figure 2 ) were found to be significantly associated with MED8 expression in HCC. In TCGA-LIHC database, using univariate analysis, we identified four immunomodulators as the significant prognostic variables for OS in HCC ( Figure 3A ). These variables were included in the stepwise multivariate Cox regression analysis, which was further used to generate a predictive model based on TGFB1, CD40LG, and TNFRSF4 ( Figure 3B ). The risk scores were calculated using the following equation: risk scores = TGFB1 × 0.155 – CD40LG × 0.8123 + TNFRSF4 × 0.295. TCGA-LIHC cohort was divided into low-risk and high-risk groups according to the optimal cut-off (0.810) obtained from the X-tile software ( Figures 3C, F ). The results of PCA and t-SNE analyses showed that low-risk and high-risk patients were distributed in two discrete groups ( Figures 3D, E ).

In TCGA-LIHC cohort, the KM curve analysis showed that patients in the high-risk group had a worse OS than those in the low-risk group ( Figure 3G ). The AUCs of the model were 0.684, 0.711, and 0.691 for 1-, 3- and 5-years OS curves, respectively, which was better than the TNM model, indicating the predictive model had a considerably good value in predicting the OS of HCC patients ( Figures 3H–J ). ROC analyses, along with the risk scores and the TNM model, increased the AUCs to 0.722, 0.756, and 0.720 for 1-, 3- and 5-years of OS, respectively. Moreover, our model showed superior accuracy in predicting long-term survival including 3- and 5-year OS relative to the previously published immune-related models ( Figures 4A –C ). There was a positive correlation between risk scores and multiple clinicopathological features ( Figures 4D –F ). The stratified survival analyses revealed that patients in the high-risk subgroup had poor OS in T1 stage ( Figure 4G ), N0 stage ( Figure 4H ), M0 stage ( Figure 4I ) and pathological stage 1 ( Figure 4J ). No significance was observed in histological grade 1, probably due to the small amount ( Figure 4K ).

Next, we determined whether this model could be used to predict the DFS. The KM curve analysis showed that the DFS for patients with high-risk scores was shorter than for those with low-risk scores ( Figure 4L ). The AUCs of the model were 0.645, 0.628, and 0.621 for 1-, 3- and 5-years DFS, respectively. ( Figure 4M ).

Furthermore, the predictive performance of the model was validated in the ICGC cohort to validate its predictive performances. The risk scores, OS, and risk-gene expression in the ICGC cohort are presented in the dot plot ( Figure 5A ). Likewise, PCA and t-SNE analyses showed that the patients were bimodally distributed ( Figures 5B, C ). The KM curve analysis indicated high-risk patients had worse OS ( Figure 5D ). The AUCs were 0.643 and 0.619 for 1- and 3-years OS, respectively ( Figures 5E, F ). Similar to TCGA-LIHC cohort, a positive association of the risk scores with the pathological stage was found ( Figure 5G ).

The nomogram, after integrating the risk scores and clinicopathological features for quantitatively predicting the probabilities of 1-, 3- and 5-years OS for individual patients, was constructed ( Figure 6A ). The AUCs were 0.718, 0.745, and 0.769 for predicting in 1-, 3- and 5-year OS, respectively ( Figure 6B ). The calibration curve exhibited a marked agreement between predicted and actual outcomes in the nomogram ( Figure 6C ).

In terms of response to the ICB therapy, chemotherapy, and targeted therapy, patients in the low-risk group showed better therapy responses relative to those in the high-risk group ( Figures 6D, E ; Figures 7A–F ), and similar results were found in the ICGC cohort ( Figures 8A–G ), which suggested that the predictive model had good power in discriminating the difference in therapeutic responses among the HCC patients. The scores for type- I IFN response, type-II IFN response, B cells, neutrophils and natural killer cells were lower in the high-risk group, while that for the macrophages showed an opposite trend ( Figure 6F ). Moreover, we performed GSEA to examine the differences in the signaling-pathways between the low-risk and high-risk groups; the top three enriched pathways based on the NES scores were displayed in each group ( Figures 7G, H ). The results showed that the signaling-pathways related to drug resistance, including the extracellular matrix (ECM) receptor interaction and focal adhesion, were significantly upregulated in the high-risk group.