The large-scale cohorts in this study included RNA sequencing data and accompanying clinicopathological data for 423 HCC patients collected from The Cancer Genome Atlas (TCGA, http://www.cancer.gov) and 10 patients with HCC from the Gumz Renal cohort in the Oncomine (http://www.oncomine.com) database. The threshold was set as follows: P-value<0.05, gene rank>10% and data type: mRNA, Cancer vs. Normal.

The Affiliated Hospital of Youjiang Medical University for Nationalities (AHYMUN, Guangxi, China) cohort consisted of 107 patients diagnosed with HCC in the Department of Hepatology, Affiliated Hospital of Youjiang Medical University for Nationalities, from June 2009 to August 2018. Clinicopathological data were collected from pathology reports or electronic medical records. Samples of HCC and normal liver tissues were collected during surgery and then processed and stored at the AHYMUN tissue bank before experiments.

IHC was performed with an anti-POLD1 antibody-N-terminal (ab226848, Abcam) at a 1:500 dilution. IHC staining was conducted in accordance with the manufacturer’s instructions as previously described (19). Based on the IHC staining intensity and density, two experienced and independent pathologists/clinicians evaluated the overall IHC score (from 0 to 12), with a score of 0 to 3 indicating negative staining and a score of 4 to 12 indicating positive staining for each tissue.

In this study, an online web tool for the retrieval of interacting genes (STRING, http://string-db.org; version 10.0) was used to establish a PPI network of differentially expressed genes (DEGs) related to POLD. Then, Spearman’s correlation analysis was performed to describe the correlation between quantitative variables without a normal distribution. The gene ontology (GO) database was used for functional enrichment analyses of biological processes (BP), molecular functions (MF) and cellular components (CC). Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome pathway enrichment databases were used to assess large-scale biologic molecular datasets. The publicly available WEB-based GEne SeT AnaLysis Toolkit (WebGestalt; http://www.webgestalt.org/; Version 6.8) was used to explore the potential function of the POLD-related gene panel (20).

To further investigate the role of POLD in HCC and pan-cancers, we analyzed the abundance and frequency of POLD mutations in HCC using cBioportal for cancer genomics (http://www.cbioportal.org/). The frequency of typical gene mutations in HCC according to differential POLD1 expression was also analyzed. Significantly elevated genes between POLD1 altered and unaltered groups were screened and identified using the Limma R package.

We grouped the HCC patients (n=107, AHYMUN cohort) into POLD1^high and POLD1^low groups using the median POLD expression. Tumor Immune Estimation Resource 2.0 (TIMER 2.0, http://timer.cistrome.org/) was used to analyze the correlation between the abundance of tumor-infiltrating immune cells and POLD expression using Spearman’s test. Then, R software was used to evaluate the interactions between immune checkpoint molecules and tumor-infiltrating lymphocytes (TILs) in human cancer samples from TCGA datasets. Spearman’s test was also used to determine the relationship between POLD1 expression and tumor purity using the ESTIMATE algorithm (21).

To determine the statistical significance of differential POLD expression between tumor and normal tissues, a Student’s t-test was performed. Kaplan–Meier curves with their 95% confidence intervals (95%CIs) and log-rank tests were applied to evaluate the significance of disease-specific survival (DSS) and overall survival (OS) benefits in separate POLD expression groups and all subgroups classified by tumor microenvironment infiltration characteristics using the Kaplan–Meier Plotter (http://kmplot.com/analysis/index). Univariate and multivariate Cox regression analyses were performed to identify the proper terms to build the nomogram. A forest plot was generated with the “forestplot” R package and used to show the P-value, hazard ratio (HR) and 95% CI of each variable. A nomogram was developed based on the results of multivariate Cox proportional hazard analysis to predict the X-year overall recurrence and calculate the risk of recurrence for individual patients. Moreover, a Sankey diagram was used to explore the relationship between pathological factors and patient survival.

All statistical analyses were performed using SPSS software (version 23.0, Inc, Chicago, IL), GraphPad Prism 8.0, R software (version 3.4.3), or online web tools. All hypothetical analyses were two-sided, and P<0.05 was considered statistically significant.

This study was conducted in three stages ( Figure S1 ). To determine the most significant genes in the DDR pathway, we first screened genes involved in base excision repair, nucleotide excision repair, mismatch repair, homologous recombination, non-homologous end-joining and Fanconi anemia pathways. We arranged these genes according to the number of pathways they are involved in and found that the POLD family participated in four pathways, indicating the important value of POLDs in DDR ( Figure 1A ). Next, the differential expression in POLD family members between >12,000 tumor and normal tissues was assessed. The expression level of POLDs was significantly higher in HCC tumor tissues compared with normal tissues ( Figure 1B ). Furthermore, POLD1 had the highest Z-score in the heatmap compared with POLD2, POLD3 and POLD4, indicating its greater significance for HCC prognosis ( Figure 1C ). Based on these results, we mainly explored the important role of POLDs in DDR and the effect of POLD expression on the development and prognosis of HCC.

To further confirm the research direction, we performed Kaplan–Meier analysis of 364 patients from the TCGA cohort. High POLD1 expression remarkably predicted a worse OS of HCC patients, whereas patients with higher expression of POLD2, POLD3 and POLD4 did not exhibit significantly shorter OS compared with lower expression patients ( Figures 2A–D ). Then, we conducted subgroup survival analyses of POLD1 expression based on hepatitis virus infection, clinicopathological staging and sorafenib use in 364 patients with HCC ( Figures 2E–H ). The results suggested that elevated POLD1 expression significantly predicted poor OS in the hepatitis virus infected and uninfected subgroup. However, for HCC patients with an advanced stage, POLD1 expression was not significantly associated with prognosis (HR=1.71, P=0.13) but predicted significantly worse outcomes for those receiving sorafenib treatment (HR=4.13, P=0.026). Additionally, high POLD1 expression was remarkably related to poor PFS in 370 HCC patients ( Figure 2I ).

Next, we explored the relationship between the expression of POLD1 and the survival of HCC patients. HCC patients with higher POLD1 expression experienced a significantly increased risk of death, and the z-score of POLD1 expression confirmed that elevated POLD1 levels were associated with higher mortality ( Figure 3A ). The Kaplan–Meier curve also demonstrated that high POLD1 expression led to worse OS in 371 patients, with a median survival of 2.8 years in the POLD^high group and 6.3 years in the POLD^low group ( Figure 3B ). The high sensitivity and specificity of the independent diagnostic and prognostic value of POLD1 expression were shown by the ROC curve ( Figure 3C ). However, the accuracy decreased over time [1-year area under the curve (AUC)=0.742, 3-year AUC=0.666 and 5-year AUC=0.610], suggesting that POLD1 expression more accurately predicted the prognosis of early-stage HCC patients compared with those with advanced stages.

Next, we assessed the Gumz Liver cohort consisting of differential RNA-seq data from HCC and normal samples. The results revealed a significantly elevated POLD1 expression level in HCC compared with normal samples (n=22; P=0.043; Figure 4A ). The Sankey diagram shows the relationships between tumor grades, stages, POLD1 expression levels and survival. Low-tumor grade and early-stage HCC patients tended to have reduced expression of POLD1, whereas high tumor grades and advanced stages were related to the increased POLD1 expression ( Figure 4B ).

Based on HCC tissues and adjacent normal tissues from TCGA (n=423) and GTEx (n=533) datasets, we first compared the POLD1 expression level between pan-cancers and normal tissues. Significant differential expression of POLD1 was commonly observed between cancer tissues and adjacent normal tissues in pan-cancers, such as GMB and STAD, and the expression of POLD1 was significantly higher in the HCC tumor group compared with the normal group ( Figures S2A, B ). Second, to explore the prognostic implications of POLD1 expression in pan-cancers, we performed Cox regression analysis ( Figures S2C, D ). POLD1 expression was closely correlated with OS in several cancers, such as PRAD (HR=9.25; P=2.7e-04) and MESO (HR=2.71; P=1.7e-06), and DSS in KIRC (HR=1.93; P=7.4e-04) and ACC (HR=3.48; P=2.6e-06).

The univariate Cox regression analysis suggested that POLD1, POLD2, POLD3 and pTNM stages were closely related to the survival of HCC patients (P<0.05). Because of the interaction between these factors, multivariate Cox regression revealed a significant effect of POLD1 expression (HR=1.63, P=0.011) and pTNM stage (HR=1.61, P<0.01) on the prognosis of HCC ( Figures 5A, B ). Together, our findings have demonstrated a consistent relationship between high POLD1 expression and advanced clinicopathological staging. Thus, we included POLD1 expression and pTNM stages in a prediction model to evaluate the prognosis of HCC using a nomogram ( Figure 5C ). Finally, the nomogram provided a graphical representation of the factors, which was used to calculate the risk of recurrence for an individual patient based on the points associated with each risk factor, and the model was more accurate in short-term survival prediction (C-index=0.682, P<0.001; Figure 5D ).

Univariate Cox regression analysis indicated that POLD1, pTNM staging and grade were closely related to the survival of HCC patients (P<0.05), and multivariate Cox regression revealed the significant effects of POLD1 expression (HR=1.54, P=0.0005) and pTNM staging (HR=1.27, P=0.016) on the prognosis of HCC ( Figures S3A, B ). These findings further confirm the relationship between high POLD1 expression and advanced clinicopathological staging in HCC prognosis. Therefore, we included POLD1 expression and pTNM stage in a nomogram to assess the prognosis of HCC patients ( Figure S3C ). The model is more accurate in short-term survival prediction ( Figure S3D ).

In addition, to investigate the underlying role of POLD1, we explored the potential value of POLD1 mutation and other related genes based on multi-omics data. First, we found that the mutation was mainly due to POLD1 missense mutation on POLBc_delta and zf-C4pol ( Figure 6A ). However, the mutation frequency was only 1.1% in the TCGA cohort. Second, we analyzed the frequency of mutations in common genes in HCC according to differential POLD1 expression in 272 cases, including TP53, TTN, CTNNB1, MUC16, ALB, PCLO, RYR2, MUC4, ABCA13, APOB and POLD1. The highest mutation frequency was found in TP53 (28%). Interestingly, the POLD1^high group exhibited significantly more TP53 mutations compared with the POLD1^low group, suggesting the importance of the TP53 expression level in clinical treatment and tumor prognosis ( Figure 6B ). Then, we explored significant DEGs between the POLD1 altered and unaltered groups and found that DNAH5, TTN, TP53, NTAN1, HDAC5, TMEM51, KIAA1211, LMBR1, MCF2L and OR5L1 were markedly upregulated, whereas TP53 expression was significantly decreased in the altered group compared with the unaltered group ( Figures 6C, D ).

To explore the mechanism of POLD1 in the DDR pathway, we constructed a PPI network and identified a gene panel of critical genes, including POLD1, POLD2, POLD3, POLD4, PCNA, MSH2, MSH6, RPA1, RPA3 and LIG1 ( Figure 7A ). Additionally, we studied the correlation between these genes. Except for POLD4, the other ten hub genes showed a close linear association with each other (Spearman’s test; Figure 7B ). To analyze the function of these genes, we conducted pathway enrichment analyses, including GO (BP, CC and MF), KEGG and Reactome ( Figures 7C, D ). We identified several critical factors, such as mismatch repair and DNA replication in BP, replication fork and a protein–DNA complex in CC, damaged DNA binding and nucleotidyltransferase activity in MF, nucleotide excision repair and base excision repair in KEGG and extension of telomeres in Reactome.

High expression levels of POLD1 were found to be closely associated with tumor purity, immune scores and stromal scores in pan-cancers ( Figure 8A ). In addition, several immune cells, including B cells, T cells, dendritic cells, macrophages and neutrophils, were significantly associated with POLD1 expression to varying degrees in different cancers, especially HCC (r²>0.4, Figure 8B ). Moreover, POLD1 expression was related to the abundance of neutrophils and CD56 natural killer cells in HCC ( Figure 8C ). Meanwhile, POLD1 levels exhibited a strong relationship with most immune checkpoint molecules, including CD274, CD80, CD86, CTLA4, PDCD1 and TCGIT ( Figure 8D ). The scatter plot shown in Figure 8 demonstrates the close relationship between POLD1 expression levels and tumor purity (HR=0.141), B cells (HR=0.468), CD8⁺ T cells (HR=0.277), CD4⁺ T cells (HR=0.358), macrophages (HR=0.397), neutrophils (HR=0.364) and dendritic cells (HR=0.438). Overall, these findings suggested that POLD1 expression was correlated with the infiltration of several immune cells and reshaped the immune-excluded microenvironment.

To validate the increased expression of POLD1 in HCC samples compared with normal liver tissues, we first collected samples from HPA and explored the prognostic implications of POLD1 expression in 107 HCC patients from the AHYMUN cohort. The nuclear expression of POLD1 was significantly higher in HCC than adjacent normal tissues based on the HPA database and AHYMUN cohorts ( Figures 9A, B ). Additionally, the results suggested that increased protein expression of POLD1 was closely associated with worse OS ( P=0.018, HR=1.697) and PFS (P=0.042, HR=1.669; Figures 9C, D ).