High expression of PARD3 predicts poor prognosis in hepatocellular carcinoma

Hepatocellular carcinoma (HCC) is one of the most commonly cancers with poor prognosis and drug response. Identifying accurate therapeutic targets would facilitate precision treatment and prolong survival for HCC. In this study, we analyzed liver hepatocellular carcinoma (LIHC) RNA sequencing (RNA-seq) data from The Cancer Genome Atlas (TCGA), and identified PARD3 as one of the most significantly differentially expressed genes (DEGs). Then, we investigated the relationship between PARD3 and outcomes of HCC, and assessed predictive capacity. Moreover, we performed functional enrichment and immune infiltration analysis to evaluate functional networks related to PARD3 in HCC and explore its role in tumor immunity. PARD3 expression levels in 371 HCC tissues were dramatically higher than those in 50 paired adjacent liver tissues (p < 0.001). High PARD3 expression was associated with poor clinicopathologic feathers, such as advanced pathologic stage (p = 0.002), vascular invasion (p = 0.012) and TP53 mutation (p = 0.009). Elevated PARD3 expression also correlated with lower overall survival (OS, HR = 2.08, 95% CI = 1.45–2.98, p < 0.001) and disease-specific survival (DSS, HR = 2.00, 95% CI = 1.27–3.16, p = 0.003). 242 up-regulated and 71 down-regulated genes showed significant association with PARD3 expression, which were involved in genomic instability, response to metal ions, and metabolisms. PARD3 is involved in diverse immune infiltration levels in HCC, especially negatively related to dendritic cells (DCs), cytotoxic cells, and plasmacytoid dendritic cells (pDCs). Altogether, PARD3 could be a potential prognostic biomarker and therapeutic target of HCC.


Scientific Reports
| (2021) 11:11078 | https://doi.org/10.1038/s41598-021-90507-w www.nature.com/scientificreports/ To screen a biomarker closely related to the formation and progression of HCC, we analyzed the differential expression of PARD3 and its clinicopathological relevance, using liver hepatocellular carcinoma (LIHC) RNA sequencing (RNA-seq) data of HCC patients from The Cancer Genome Atlas (TCGA). Then, we put PARD3 into a prognosis analysis in order to evaluate predictive capacity. Moreover, we performed comparative transcriptome analysis, functional enrichment analysis and correlation analysis between PARD3 and immune cell infiltration to evaluate functional networks related to PARD3 in HCC and explore its role in tumor immunity.
Correlation of PARD3 overexpression with poor clinicopathologic features. To investigate overexpression of PARD3 and its clinicopathological relevance, 371 HCC samples with detailed patient information (retrieved from TCGA in June 2020) were divided into two groups by the median value of PARD3 expression. As shown in Table 1, high PARD3 expression was associated with advanced T stage, pathologic stage, residual tumor, histologic grade, vascular invasion and higher alpha fetoprotein (AFP). Otherwise, high PARD3 expression group also carried more TP53 mutation (Mut) than low PARD3 expression group. Whereas, the distributions of other clinicopathologic features showed no difference between high and low PARD3 expression group. Univariate logistic regression further confirmed the association between high PARD3 expression and poor clinicopathologic characteristics in HCC patients ( Fig. 2A,B). In addition, the area under receiver operation characteristic (ROC) curve (AUC, AUC = 0.835, 95% CI = 0.792-0.877) indicated that PARD3 had a good diagnostic power, and was expected to be a potential biomarker for HCC (Fig. 2C).
Correlation of PARD3 overexpression with adverse Outcomes in HCC. As PARD3 overexpression was correlated with poor clinicopathologic features, we then explored the prognostic value of PARD3 in HCC.
Based on multivariate COX regression, we constructed a nomogram to provide a visual, intuitive and quantitative description of those risk factors and their weights in HCC, and predict the probability of 1-year, 3-year and 5-year survival. Potential covariates included adjacent hepatic tissue inflammation (39 points), T stage (75.5 points), age (67 points), tumor status (100 points) and PARD3 expression (90.25 points), and higher total points indicated worse prognosis (Fig. 3C).
Then, we made calibration curve of the nomogram, and calculated concordance index (C-index) to assess the predictive ability of PARD3 as a biomarker for HCC. As shown in Fig. 3D, the bias-corrected line in the calibration curve was close to the ideal line, and the C-index was 0.702 (95% CI = 0.668-0.736). In all, the nomogram was available to predict the prognosis of HCC patients, and PARD3 exhibited stable predictive ability.
The result of CCs indicated that PARD3 and its associated DEGs were primarily located in lipoprotein particles, synaptic membranes and neuron projection membranes (Fig. 4B). In terms of BPs and MFs, these genes mainly participated in detoxification of copper ion (Cu), stress response to metal ions, and zinc ion (Zn) homeostasis, which associated with receptor ligand activity, metal ion transmembrane transporter activity, DNA-binding transcription activator activity, hormone activity, carbohydrate binding, and arachidonic acid oxygenase activity (Fig. 4C,D).
Signaling pathways, working as functional units of gene groups, play an important role in cell biological effects. Thus, we identified significantly enriched signaling pathways between low and high PARD3 expression groups by Gene set enrichment analysis (GSEA). According to normalized enrichment scores (NSEs), 431 pathways were significantly associated with high expression of PARD3, including cell cycle and mitosis, DNA double strand break repair, cell motility, Rho, MAPK, TP53, response to metal ions, and pathways related to metabolism ( Table 2). In addition, we made a protein-protein interaction (PPI) to highlight the most important protein functional groups interacting with each other. The two most crucial MCODE subnetworks were SAA1-related cluster and CYP-related cluster, both of which were involved in metabolism and homeostasis (Fig. 5). The above result suggested a negative impact of PARD3 on tumorigenesis and progression of HCC. www.nature.com/scientificreports/ of 24 types of tumor-infiltrating immune cells (TIICs), in order to evaluate the association between PARD3 and immune infiltration levels in HCC. As illustrated in Fig. 6, PARD3 was involved in infiltration of T helper cells, Th2 cells, and T central memory (Tcm); but negatively related to infiltration of dendritic cells (DCs), cytotoxic cells, plasmacytoid dendritic cells (pDCs), neutrophils, immature DCs (iDCs), and regulatory T cells (Treg). Furthermore, we replicated immune infiltration analysis using another tumor-immune system interaction database (TISIDB), and obtained consistent results (Fig. S2).

Discussion
PARD3 plays a crucial role in establishment and maintenance of epithelial cell polarity 23 . So far, at least five PARD3 variants have been identified in human liver cDNA library 32 . PARD3 largely engages in cancer cell proliferation, apoptosis, invasion, migration and epithelial-mesenchymal transition (EMT) 18,24,33,34 . The modulation of PARD3 in tumorigenesis and progression among different cancers seems to be controversial. For instance, PARD3 acts as a tumor suppressor in lung, bladder, breast, cervical, esophageal and pancreatic cancers and malignant melanoma 15,22,24,27,[35][36][37][38] , but it is found to be activated in ovarian cancer and clear-cell renal carcinoma 18,25,39 . In skin cancers, PARD3 shows dual effects depending on the tumor type 16 . In terms of HCC, a study reported the association between overexpression of PARD3 and extrahepatic metastasis, and suggested one of its possible mechanisms. However, the specific role and detailed mechanism of PARD3 in HCC has not been fully elucidated 17 . Hence, we performed the bioinformatic analysis using several independent databases Based on TCGA database, PARD3 expressed differentially in 24 types of cancers (Fig. 1A). Thereinto, the expression levels in two types of cancers were inconsistent with previous studies, including lung cancer and pancreatic cancer 22,27 , which may derive in part from data collection approaches and patients' biological properties. More specifically, our study provided direct evidences that PARD3 was an independent risk factor for HCC development. Firstly, PARD3 expression was dramatically higher in HCC than in normal liver tissues (Fig. 1B-I). Secondly, high PARD3 expression group contained more patients with advanced pathologic stages, vascular invasion, and TP53 Mut, suggesting that overexpression of PARD3 was significantly associated with poor clinicopathologic features with a good predictive power (AUC = 0.835, Fig. 2 and Table 1). Thirdly, elevated PARD3 expression led to shorter OS and DSS in both whites and Asians regardless of gender, age and Child-Pugh grade (Fig. 3A,B). All these results proved PARD3 as a potential prognostic biomarker of HCC.
Thus, we further explored the possible mechanism by which high PARD3 expression worsens the outcomes of HCC patients. As an adverse prognostic indicator of HCC, PARD3 was involved in many pivotal mechanisms in cancer, including cell cycle 40 , DNA damage and repair 41 , and cell motility 42 (Table 2). It is well known that genomic instability and mutagenesis, which caused by erroneous DNA repair, are closely correlated with poor prognosis and drug resistance in HCC 41 . TP53 is universally recognized as a hub gene in responding to DNA damage and guarding the genome, and its mutation is observed in about half of all solid tumors, including HCC 43,44 . Furthermore, p38MAPK is able to control p53 activation via direct phosphorylation. Based on our result of enrichment analysis that PARD3 was associated with MAPK pathway and TP53 regulation (Table 2), we hypothesized PARD3 may affect the formation and progression of HCC by regulating TP53 via MAPK pathway. www.nature.com/scientificreports/ However, the hypothesis requires further investigation. Rho Family GTPases, which is closely interact with PARD3, were widely reported to regulate cell cycle and cell motility across human cancer of different origins 42 . PARD3 directly activates Rac1, promoting proliferation and motility of cancer cells, and leads to tumorigenesis, angiogenesis, invasion and metastasis [45][46][47] . Likewise, Par complex also links Rho small GTPases to regulate asymmetrical cell division and cell polarization 48 , which manipulate EMT and mesenchymal-epithelial transition (MET) 49 . Our findings that PARD3 was significantly implicated in Rho pathway also provided evidences to confirm this theory (Table 2). In addition, some lncRNAs showed significant correlations with PARD3 expression, such as FAM83A-AS1, which is involved in HCC (Table 2) 50 .
Remarkably, we found that PARD3 and its associated DEGs mainly participated in cellular response to metal ions (Fig. 4C,D, and Table 2). Previous studies have provided a possible relationship between metal ion homeostasis and vascular invasion in HCC, which may be mediated by p53 [51][52][53] . Disturbance in Cu and Zn homeostasis has been reported as a significant factor associated with tumor proliferation, angiogenesis and invasion in HCC, and furthermore, cellular response to Cu and Zn is probably involved in mitochondrial accumulation and stability of p53, so as to influence proliferation and apoptosis of hepatoma cells [54][55][56] .
In the past decades, reprogramming of energy metabolism was added into the list of cancer hallmarks 29 . Interestingly, we found many co-expression genes and pathways related to PARD3 were involved in deregulation of metabolisms, covering types of metabolic processes like fat acids, amino acids and pyrimidines ( Table 2). In order to sustain prodigious proliferation, tumors exert a specialized metabolism that differs from normal tissues. During the period, tumors recruit abundant nucleotides to maintain unlimited replicative potential, and uptake more nutrients to support unchecked cell growth. In particular, alterations in metabolism fatty acid and glycine, serine and threonine have been investigated as a promoter of HCC initiation and progression 57,58 . Moreover, SLC22A1, a DEG inversely related to PARD3 (Fig. 4E), is a key regulator of metabolism, which is extensively considered as a suppressor of HCC development [59][60][61] .
Besides, PPI enrichment analysis screened SAA1-related cluster and CYP-related cluster as the two most crucial MCODE subnetworks, both of which were involved in metabolism and homeostasis (Fig. 5). Recent study showed that downregulated SAA1 was closely associated with progression of HCC and low anti-tumor immune infiltrating 62 ; and CYP families might impact HCC cell viability via modulating biotransformation 63 . The results provided supporting evidence that PARD3 might promote HCC via regulating metabolism and homeostasis.
Recently, cellular metabolism has emerged as a determinant of the viability and function of both tumor cells and immune cells. Meanwhile, tumor metabolism is reported as an immune checkpoint 58,64 . As discussed above, PARD3 was linked to some important metabolic processes, and meanwhile, several enriched pathways were also associated with immune response (Table S4). Thus, we hypothesized that there may be an association between PARD3 and immune infiltration. As expected, PARD3 is correlated with diverse immune infiltration levels in HCC, especially DCs, cytotoxic cells and pDCs (Fig. 6). DCs are a heterogeneous population of professional antigen-presenting cells central to the induction and maintenance of adaptive immunity within tumor microenvironment 58,65 . In particular, two subsets of DCs exert the most potent antitumor functions, including conventional DCs type 1 (cDC1s) that stimulate T cell proliferation, and pDCs that produce interferon-α (IFNα) 65,66 . cDC1s not only take up and cross-present tumor antigens via major histocompatibility complex (MHC) class I to activate naive CD8 + T cells; but also support the cytotoxicity of CD8 + T cells by secreting large amounts of interleukin-12 (IL-12). Then, activated cytotoxic CD8 + T cells migrate to tumors and kill them 65 . pCDs play two opposite roles in tumor immunity depending on their subsets via inducing Treg or activating cytotoxic T cells respectively 65,66 . Based on our result that PARD3 negatively correlated with DCs and cytotoxic T cells, we speculate that immune infiltration related to PARD3 may contribute to the unfavorable outcomes for HCC, yet the specific regulation mechanism needed to further elucidate.
In summary, our study reveals that overexpression of PARD3 correlates with poor clinicopathologic features and adverse outcomes in HCC. Moreover, the crosstalk of cellular response of metal ion, metabolism and immune infiltration within tumor microenvironment may partly explain the function of PARD3 in HCC development. However, bioinformatic analysis based on TCGA also has some limitations. First, the sample sizes of blacks and stage IV in LIHC may be too small to show a significant difference between groups. Additionally, transcriptome sequencing cannot directly reflect the protein activity and expression level. Therefore, our results should be verified by further research using sufficient HCC clinical samples, and detailed mechanisms of PARD3 need investigating more intensively. Despite the limitations, our findings provide multilevel evidence for the value of PARD3 as a potential prognostic biomarker and therapeutic target of HCC.

Materials and methods
RNA sequencing data and processing. RNA-Seq data (Workflow Type: HTSeq-FPKM) and corresponding clinical information were retrieved from TCGA-LIHC project, among which 371 HCC patients with complete survival information were retained. Then, level 3 HTseq-FPKM data were transformed to transcripts per million reads (TPM) for further analysis. Unavailable or unknown clinical data were treated as missing values 67 . RNA-Seq data of multiple cancer types were downloaded from the online database UCSC Xena (https:// xenab rowser. net/ datap ages/), and analyzed using Toil 68 . This study complied with the publication guidelines provided by TCGA.
Differentially expressed gene analysis. DESeq2 package was used to identify DEGs 69 . The cut-off value of PARD3 expression was determined by its median value, and the thresholds were defined as |log fold change (log FC)| > 2 and adjusted p-value < 0.05. www.nature.com/scientificreports/ The differential expression of PARD3 was simultaneously validated using Gene Expression Omnibus (GEO) database (https:// www. ncbi. nlm. nih. gov/ geo), including three independent cohorts (GSE14520, GSE76427 and GSE121248) 70-72 . Functional enrichment analysis. Metascape analysis. Metascape (http:// metas cape. org) was used as a gene list analysis tool to conduct GO enrichment analysis of DEGs 73 , including BPs, MFs and CCs. P-value < 0.05, minimum count > 3 and enrichment factor > 1.5 were considered to be significant. The Cytoscape plugin MCODE was used to screen crucial clustering subnetworks in PPI network.
Gene set enrichment analysis. GSEA was used for Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment, which was performed with 1000 permutations for each analysis using curated gene sets (C2. cp.v7.0.symbols.gmt) as the reference gene set 74 . Visualization and statistics were carried out by R package clusterProfiler 75 . Adjusted p-value < 0.05, false discovery rate (FDR) q-value < 0.25 and |normalized enrichment score (NES)| > 1 were considered to be significant.
Immune infiltration analysis. The relative abundance of each immunocyte type was described with EC in single-sample Gene Set Enrichment Analysis (ssGSEA). ECs for 24 types of TIICs were quantified using GSVA package in R as reported previously 76 .
Prognostic model generation and statistical analysis. All statistics were performed using R (v.3.6.2).
The PARD3 expression levels between tumor and normal tissues in paired or non-paired samples were compared using Wilcoxon signed-rank test and Wilcoxon rank sum test, respectively. The discrimination ability of PARD3 in HCC was evaluated using the AUC in ROC 78 . The correlation between PARD3 expression and screened DEGs were analyzed using Spearman's correlation. The correlation between PARD3 expression and immune infiltration were analyzed using Spearman's correlation, while ECs of the two groups with different expression level were compared using Wilcoxon rank test. The relationship between clinicopathological features and PARD3 expression in HCC patients was assessed using Kruskal-Wallis test, Wilcoxon rank test or Spearman's correlation, while prognostic relevance of the two groups with different expression level were compared using Pearson χ 2 test, Fisher exact test or univariate logistic regression. A survival curve was plotted using Kaplan-Meier method, and analyzed by Cox regression. Baseline variables with a p-value < 0.1 on univariate analysis were included in multivariate Cox regression model 79,80 . A nomogram was generated to predict the prognosis of HCC based on the result of multivariate Cox regression analysis, including significant clinical characteristics and PARD3 expression. C-index was used to validate the predictive power of the model 81 . Statistical results were displayed with p-value, and hazard ratio (HR) at a 95% confidence interval (95% CI). p-value < 0.05 were considered to be statistically significant.

Data availability
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.