Elevated DNA Polymerase Delta 1 Expression Correlates With Tumor Progression and Immunosuppressive Tumor Microenvironment in Hepatocellular Carcinoma

Background and Objective Hepatocellular carcinoma (HCC) is one of the most common cancers worldwide, and the DNA polymerase delta (POLD) family is significantly related to cancer prognosis. This study aimed to explore the significance of the POLD family in HCC via the DNA damage repair (DDR) pathway. Methods Data mining was conducted using bioinformatics methods. RNA sequencing and clinicopathological data were collected from The Cancer Genome Atlas, GTEx database and the Gumz Renal cohort. Statistical analyses were also performed in cancer samples (n>12,000) and the Affiliated Hospital of Youjiang Medical University for Nationalities (AHYMUN, n=107) cohort. Results The POLD family (POLD1–4) was identified as the most important functional component of the DDR pathway. Based on the analysis of independent cohorts, we found significantly elevated POLD expression in HCC compared with normal tissues. Second, we investigated the prognostic implication of elevated POLD1 expression in HCC and pan-cancers, revealing that increased POLD1 levels were correlated to worse prognoses for HCC patients. Additionally, we identified 11 hub proteins interacting closely with POLD proteins in base excision repair, protein-DNA complex and mismatch repair signaling pathways. Moreover, POLD1 mutation functioned as an independent biomarker to predict the benefit of targeted treatment. Importantly, POLD1 expression was associated with immune checkpoint molecules, including CD274, CD80, CD86, CTLA4, PDCD1 and TCGIT, and facilitated an immune-excluded tumor microenvironment. Additionally, we confirmed that elevated POLD1 expression was closely correlated with the aggressive progression and poor prognosis of HCC in the real-world AHYMUN cohort. Conclusion We identified a significant association between elevated POLD1 expression and poor patient survival and immune-excluded tumor microenvironment of HCC. Together, these findings indicate that POLD1 provides a valuable biomarker to guide the molecular diagnosis and development of novel targeted therapeutic strategies for HCC patients.


INTRODUCTION
Hepatocellular carcinoma (HCC) is the most common type of primary liver cancer, with over half a million new cases diagnosed worldwide each year. HCC is the most malignant form of liver cancer, which derives from the aggressive proliferation of liver epithelial cells (1). As the sixth leading cause of cancer-related mortality globally, HCC accounts for 4.7% of all cancer deaths. In 2020, over 910,000 people were diagnosed with HCC, and its incidence is continuously rising (2). In China, the incidence of HCC is particularly high, accounting for 55% of the total number of HCC patients worldwide.
Several molecular pathways are implicated in HCC carcinogenesis, including vascular endothelial growth factor receptor (VEGFR), TP53 and Akt/mTOR pathways (3). Among several pathogenic factors, cirrhosis underlies most cases of HCC, and hepatitis C virus (HCV) infection also causes cirrhosis and HCC. However, the eradication of HCVrelated cirrhosis does not prevent the development of HCC (4,5). Thus, more effective HCC treatments are needed.
Historically, serum alpha-fetoprotein (AFP) and diagnostic imaging were the primary diagnostic methods for HCC. However, the poor prognosis of patients owing to the late diagnosis of HCC is unacceptable, and AFP levels are not significantly elevated in most small and early-stage HCCs (6). Recently, several studies have focused on identifying promising biomarkers for early detection (7,8). Currently, targeted drugs, such as sorafenib, are used as a first-line treatment strategy for the personalized treatment of advanced HCC (9). However, exploring unknown targets and pathways is important to achieve more significant survival benefits for advanced HCC patients who cannot be cured through surgery.
DNA damage repair (DDR) pathways consist of multiple interconnected cellular signaling networks that are activated in response to DNA damage. DDR is correlated with genomic instability, tumor mutational burden in HCC and immune cell function (10). Several studies have shown that the development of cancer is mainly driven by defects in DDR. Therefore, inhibitors targeting kinases involved in DDR, such as AZ20 and M3814, have been investigated for the treatment of cancers in clinical trials (11).
The DNA polymerase delta (POLD) family consisting of POLD1, POLD2, POLD3 and POLD4 is an important mediator of DNA repair during chromosome replication, which is of great significance to the maintenance of DNA structure stability (12). POLD1 encodes the catalytic subunit of DNA polymerase delta, which is involved in both polymerase and 3' to 5' exonuclease activities. POLD1 plays a crucial role in DNA replication and DNA double-strand break repair, together with the accessory proteins POLD2, POLD3 and POLD4 for full activity (13,14). POLD2 encodes the 50-kDa catalytic subunit of DNA polymerase delta, assisting POLD1 in the process of DDR.
In addition, POLD3 and POLD4 may catalyze the repair of damaged replication forks via break-induced replication (15). Recent studies have shown that the loss or reduced expression of POLD3 predicts the occurrence and progression of colorectal cancer and other malignant tumors (16,17). The decreased expression of POLD3 is associated with the pathogenesis of HCC by enhancing the proliferation and invasion abilities of tumor cells (18).
However, the relationship between the expression of POLD family members and the prognosis of HCC remains poorly understood. Based on the role of the POLD family in the maintenance and repair of DNA, we hypothesized that POLD expression significantly affects the prognosis of HCC via the DDR pathway. Our study aims to identify a biomarker with significant potential to improve the diagnosis and prognosis of HCC.

Patients and Tissue Samples From Online Databases and Real-World Cohorts
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.

Immunohistochemistry (IHC) Staining Analysis
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.

Establishment of a Protein-Protein Interaction (PPI) Network and Functional Annotations of POLD-Related PPI Networks
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).

Abundance and Frequency of POLD Mutations in HCC
To further investigate the role of POLD in HCC and pancancers, 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.

Analysis of Immune Infiltration in the Tumor Microenvironment
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).

Statistics Analysis
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.

Identification of Critical DDR Pathway Genes and Their Prognostic Implications in HCC
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.

Prognostic Role of POLDs in 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).

POLD1 Significantly Predicts Survival and Aggressive Clinicopathological Parameters for HCC Patients
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 The expression levels of POLDs were significantly higher in HCC tumor tissues compared with normal tissues in more than 12,000 samples. *p < 0.05, **p < 0.01, ***p < 0.001. (C) Heatmap analysis indicated that POLD1 has the most significant influence on HCC prognosis compared with POLD2, POLD3 and POLD4.

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).

Differential Expression of POLD1 in Pan-Cancers and Its Prognostic Value
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).

Construction of a Prediction Model Using Cox Regression and Nomogram Analysis
The univariate Cox regression analysis suggested that POLD1, POLD2, POLD3 and pTNM stages were closely related to the   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 shortterm survival prediction ( Figure S3D).

Landscape of POLD1 Mutation and Related Genes in HCC
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).

PPI Network Establishment and Functional Enrichment Analysis
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.

POLD1 Expression Predicts the Abundance of Immune Cells and Immune Checkpoint Molecules in the HCC Microenvironment
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 2 >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  Overall, these findings suggested that POLD1 expression was correlated with the infiltration of several immune cells and reshaped the immune-excluded microenvironment.

Validation of Differential POLD1 Expression and Its Prognostic Value in the AHYMUN Cohort
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).

DISCUSSION
HCC is the most common primary liver cancer and has a poor prognosis. Over the past few decades, the incidences of liver cancer and liver cancer-related deaths have increased (4,22). However, the treatment options for advanced liver cancer are very limited, and strategies for advanced personalized treatment of HCC are lacking (23). Therefore, better understanding of the mechanisms underlying HCC is critical because these findings may help identify novel treatments for HCC patients. In this study, we examined the potential prognostic value of POLD in HCC and evaluated its role in the tumor immune microenvironment. POLD1 interacts with POLD2, POLD3, and POLD4 to form the POLd holoenzyme together with replication factor C and proliferating cell nuclear antigen. POLD2-4 are smaller subunits, and POLD1 plays a major role in the biochemical activity of the polymerase (24). Consistent with these data, we found that POLD1 showed a more significant influence on tumor prognosis than the other three POLD proteins and therefore we focused on POLD1 for the subsequent analysis.
POLd has important proofreading capabilities conferred by exonuclease activity and is also involved in repairing DNA damage, including nucleotide excision repair, double-strand break repair, base excision repair, and mismatch repair (25,26). Multiple studies have linked germline and sporadic mutations in POLD1 and other subunits of POLd with human pathologies (27,28). For example, mutations in POLd in mice and humans lead to genome instability, mutant phenotypes, and tumorigenesis. In addition, mutations in the proofreading domain of POLD1 have been identified as the root cause of some hereditary cancers and these mutations may affect treatment management. Recent studies have shown that lossof-function alterations in DDR genes are associated with human tumorigenesis (29). Germline mutations in POLE and POLD1 have been shown to predispose patients to multiple colorectal polyps and a wide range of neoplasms (30). A previous study showed that at least 1 in 92 primary liver cancer patients had DDR gene mutation (31). The landscape of DDR mutations and their association with genetic and clinicopathologic features suggest that PLC patients with altered DDR genes may be rational candidates for precision treatment (32). Therefore, in addition to the increased expression level of POLD1, POLD1 and POLE mutations may function as independent biomarkers to predict the benefit of targeted treatment.
The expression levels of POLD family members in tumor and normal tissues were explored. According to previous research, the increased expression of POLD3 indicated a poor prognosis of HCC patients, although not as significantly as POLD1 (33). POLD3 plays a specialized role in the repair of damaged replication forks, indicating that POLD3 activity may be particularly relevant for cancer cells enduring high levels of DNA replication stress (12). The cellular depletion of POLD1 or POLD3 resulted in differential genome instability manifested by DNA double-stranded breaks (17,34).
Human cancers can be divided into three types according to the anti-tumor immune response status, or the immune phenotype: inflamed, immune-excluded, and immune desert status. Inflamed cancers generally refer to tumors with high PD-L1 expression in cancer cells and more immune cells and tumor infiltrating lymphocytes in the tumor; these tumors are sensitive to immune checkpoint inhibitors. Immune-excluded tumors are tumors in which the stroma shows a large number of T cells, but the T cells cannot penetrate the stroma and infiltrate the tumor because of the strong inhibitory microenvironment. These tumors often do not respond well to immune checkpoint inhibitors. Immune desert tumors lack the infiltration of T cells and immune cells even in the interstitium, which is a described as an "immune desert". In our analysis, we found that POLD1 is closely related to the expression of B cells, CD8+ T cells, CD4+ T cells and macrophages. These results suggest that POLD1 may play an important role in the immune-excluded tumor microenvironment. We speculate that this may be because POLD1 is involved in DNA cleavage repair, which affects the tumor microenvironment.
There are some advantages of this study. First, this research contained independent HCC cohorts, including the TCGA database (n=423), GTEx database and Gumz Renal cohort (n=10), and the real-world AHYMUN cohort (n=107). Second, we demonstrated the value of significantly elevated POLD expression for HCC prognosis and identified POLD1 as the most valuable gene for further analysis. Third, at both the mRNA and protein levels, we validated the association between POLD1 expression and HCC prognosis. Moreover,~12,000 tumor samples from TCGA database were collected, and the effect of POLD1 expression on prognosis was verified in pan-cancers. Finally, functional enrichment analysis was performed. The role of POLD1 in the infiltration of immune cells in the tumor microenvironment was demonstrated, which may guide cancer treatment and targeted drug development. This study has several limitations. First, there is a high degree of heterogeneity between the patient groups in this study. Therefore, in the next study, we will select patients from multiple regions for multi-center research. Second, this was a retrospective study. Finally, we will further study the demographic, clinical and pathological details of the population.
Our study has demonstrated a link between elevated POLD1 expression and patient survival and the tumor microenvironment in HCC. In the next study, we will explore the role and mechanism of POLD in HCC progression through cytological tests and animal experiments. We will also analyze the potential mechanisms of POLD, POLd, and DNA cleavage repair in HCC through high-throughput sequencing and other methods. These findings will help provide new insights into the pathogenesis of HCC and new ideas for the personalized treatment of HCC patients.

CONCLUSION
This study first investigated the molecular and clinical role of the POLD family and revealed the significant relationship between elevated POLD1 expression and the poor survival and immuneexcluded tumor microenvironment of HCC patients. Together, these findings support the use of POLD1 as a biomarker to guide the molecular diagnosis and development of novel targeted therapeutic strategies for HCC patients.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

ETHICS STATEMENT
The study design and experimental procedures were performed in accordance with the Declaration of Helsinki. The study was approved by the ethics committee of the Affiliated Hospital of Youjiang Medical College for Nationalities (Baise, Guangxi Province, China). Written informed consent was obtained from all participants.

AUTHOR CONTRIBUTIONS
WL and SZ carried out the molecular genetic studies, participated in the sequence alignment and drafted the manuscript. JW and SC carried out the immunoassays. HD, CW and BH participated in the sequence alignment. XS and QG participated in the design of the study and performed the statistical analyses. HT, SZ, YZ, WL, and JW conceived the study, participated in the study design and coordination and helped to draft the manuscript. All authors contributed to the article and approved the submitted version.

SUPPLEMENTARY MATERIAL
The The prediction is more accurate in short-term survival prediction.