Overexpressed PAQR4 predicts poor overall survival and construction of a prognostic nomogram based on PAQR family for hepatocellular carcinoma

1.   Introduction

In 2020, the global cancer report showed 905,677 new cases and 830,180 deaths of liver cancer; additionally, the sixth highest incidence and fourth highest mortality rates were observed for the common malignant tumors of the liver ^[1]. According to histology, primary liver cancer is divided into hepatocellular carcinoma (HCC) (approximately 75–85%), intrahepatic cholangiocarcinoma (ICC) (approximately 10–15%), and other rare types ^[2]. Although surgical resection is considered as the optimal treatment for patients at an early stage, most patients are diagnosed at an advanced stage ^[3]. However, the recurrence and metastasis of HCC cannot be resolved. Currently, the overall prognosis of patients with liver cancer is poor, with a 5-year overall survival rate of approximately 18% ^[4]. In recent years, as a new strategy, immunotherapy has played a pivotal role in advanced HCC, among which immune checkpoint inhibitors (ICIs) have been nearly the front-line treatment in advanced HCC and significantly improved clinical outcomes of advanced patients ^[5,6]. Currently, ICIs as monotherapy have shown controversial results in HCC, while combination therapy between ICIs has achieved more surprising therapeutic effects, such as lenvatinib plus pembrolizumab and atezolizumab plus bevacizumab ^[7,8]. Therefore, substantial efforts are still required to identify the novel biomarkers that contribute to the improved effectiveness of the diagnosis, treatment, and prognosis of HCC.

In the human genome, there are 11 members of the PAQR family (PAQR1 to PAQR11), including PAQR4, which are mainly involved in regulating metabolism and carcinogenesis ^[9,10]. PAQR1-4 include adiponectin-related receptor subtype, in which AdipoR1 (PAQR1) and AdipoR2 (PAQR2) regulate fatty acid oxidation and glucose uptake ^[11]. PAQR3 plays a role in cell cycle of tumor cell proliferation and apoptosis ^[12]. In humans, PAQR5-9 include membrane progesterone receptors (mPRs) containing the five sub-types, PAQR5 (mPRγ), PAQR6 (mPRδ), PAQR7 (mPRα), PAQR8 (mPRβ), and PAQR9 (mPRε), which regulate the cell cycle processes and tumor malignant biological behaviors ^[10,13]. MMD2 (PAQR10) and MMD (PAQR11) are hemolysin-related receptor subtypes that mediate the Ras signaling pathway in the Golgi apparatus ^[14]. An increasing number of studies have revealed that the upregulation of PAQR4 is related to the occurrence, progression, and adverse prognosis of multiple malignancies ^[15,16,17]. Nevertheless, the carcinogenic effect and prognosis of PAQR4 in HCC have not yet been investigated.

In the present study, we analyzed the expression level of PAQR4, its correlation with clinical characteristics, methylation level, prognostic value, immune implication, and the potential pathways of PAQR4 in HCC. At the same time, the prognostic significance of PAQR family genes in HCC was analyzed using Cox regression; additionally, a clinical prognosis nomogram model was constructed by combining the prognostic genes and clinicopathological parameters.

2.   Materials and methods

2.1.   Data collection

A flow-process diagram of the study design is shown in Figure 1. HCC datasets, including RNA-seq data (including 374 tumor samples and 50 normal samples), clinicopathological data, and DNA methylation data, were collected from the TCGA database (https://portal.gdc.cancer.gov/; accessed April 20, 2021). The HCC datasets, which include RNA-seq data (including 232 tumor samples and 202 normal samples) and clinical data were collected from the ICGC database (https://dcc.icgc.org/; accessed April 20, 2021). GSE14520 (including 225 tumor samples and 220 normal samples) and GSE36376 (including 240 tumor samples and 193 normal samples) datasets of RNA-seq data were collected from the GEO database (https://www.ncbi.nlm.nih.gov/geo/; accessed April 21, 2021). The gene expression matrices were processed by the "limma" package of R software ^[18].

Figure 1.  The work flow chart of this study. Abbreviations: TCGA, The Cancer Genome Atlas; ICGC, International Cancer Genome Consortium; GEO, Gene Expression Omnibus; GSEA, Gene Set Enrichment Analysis; ROC, Receiver Operator Characteristic.

DownLoad: Full-Size Img PowerPoint

2.2.   Association between PAQR4 and clinicopathological features

The difference in the expression of PAQR4 between HCC tissues and normal liver tissues was analyzed using multiple databases. Clinical HCC samples with incomplete survival time and survival status were excluded, and 367 samples and 232 samples were retained from the TCGA and ICGC datasets, respectively. We divided the HCC patients into high PAQR4 expression group and low PAQR4 expression group according to the median PAQR4 expression value. Subsequently, Kaplan-Meier curves were plotted using the ''survminer'' package and ''survival'' package of R software. In addition, the correlation between the clinicopathological factors of HCC patients and PAQR4 expression was analyzed using the TCGA dataset.

2.3.   Copy number variation, DNA methylation and regulatory miRNA of PAQR4

The relationship between PAQR4 expression and copy number variation (CNV) was analyzed using the cBioPortal database (http://www.cbioportal.org/) ^[19]. The correlation between PAQR4 methylation and PAQR4 expression was analyzed using the Pearson correlation coefficient. The miRNAs that regulated PAQR4 were analyzed using the TargetScanv7.2 database (http://www.targetscan.org/vert_72/) and miRWalk database (http://mirwalk.umm.uni-heidelberg.de/) ^[20,21]. The downregulated miRNAs in HCC were screened from the TCGA dataset. The miRNAs shared among the above three databases were considered as the potential regulators of PAQR4 in HCC.

2.4.   Immune cells and immune checkpoint molecules

TIMER (https://cistrome.shinyapps.io/timer/) ^[22] was utilized to analyze the association between the tumor-infiltrating immune cells (TIICs) and immune checkpoint molecules with PAQR4.

2.5.   Gene set enrichment analysis

We explored the PAQR4-related signaling pathways in HCC using GSEA (version 4.1.0). The "c2.cp.kegg.v7.2.symbols.gmt" was identified as the reference gene set. Gene set permutation was performed 1000 times for each analysis. An enrichment plot was drawn using the "plyr" package, "ggplot2" package, "grid" package and "gridextra" package. Significant enrichment pathways were identified with P < 0.05, and false discovery rate (FDR) < 0.05 ^[23].

2.6.   Cell culture

Five hepatoma cell lines (HEPG2, LM3, SKHEP1, MHCC97H and HUH7) and normal hepatocyte LO2 were obtained from the Beijing Cell Resource Center. The cells were cultured in RPMI1640 medium containing 10% fetal bovine serum, 5% CO₂ and 37 ℃ in a constant temperature and humidity incubator.

2.7.   qRT-PCR

Six pairs of clinical HCC tissues and adjacent tissues were obtained from patients with primary liver cancer; additionally, all patients underwent radical hepatectomy in the department of general surgery. This study was approved by the institutional research ethics board of Lanzhou university second hospital (protocol number: 2019A-321) and conformed to the ethical ethics of the 1975 Declaration of Helsinki. Total RNA was extracted from cultured cell lines and frozen tissue specimens using TRIzol reagent (Invitrogen) according to the manufacturer's instructions; the purity and concentration of the extracted RNA were detected using a spectrophotometer; and the RT reagent kit (Accurate Biology, China) was used to reverse transcribe RNA into cDNA. We designed the primer sequence of PAQR4 (Forward: 5'-TGCCCGCCTACTGGTATTTG-3' and Reverse: 5'-GCTCAGCAGGTGCATGATCT-3') and GAPDH (Forward: 5'-GGAGCGAGATCCCTCCAAAAT-3' and Reverse: 5'-GGCTGTTGTCATACTTCTCATGG-3') using the PrimerBank database (https://pga.mgh.harvard.edu/primerbank/) ^[24]. Further, with the aid of NCBI BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) ^[25] it was shown that the sequence only targets on PAQR4. Amplification was performed using a CFX96 PCR system (Bio-Rad, United States). The reaction conditions were as follows: predenaturation at 95 ℃ for 30 s, followed by denaturation at 95 ℃ for 5 s, and a final annealing and extension at 60 ℃ for 30 s were carried out, for a total of 40 cycles. Relative quantitative data were calculated with the 2^-ΔΔCT method.

2.8.   The PAQR family

The PAQR family has 11 members in humans, including PAQR1-11. Protein-protein interaction (PPI) of the PAQR family was analyzed using STRING network tools (https://string-db.org/) ^[26]. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of PAQR family genes were performed using the David database (https://david.ncifcrf.gov/) ^[27]. GO includes biological process (BP), cellular component (CC) and molecular function (MF). In addition, correlations between PAQR family genes were analyzed by the Pearson correlation coefficient. The gene expression heatmap of the PAQR family in HCC was plotted using the "pheatmap" package.

2.9.   Prognostic model construction

We identified the prognostic gene signatures of the PAQR family in HCC using Cox regression analysis and established a credible prognostic model. Risk score = ∑(coefficient_i × expression value of gene_i). We divided the patients with HCC into high-risk and low-risk groups based on median value of the risk score. The survival differences between the two groups were analyzed using the Kaplan-Meier curve. We used the ''plotROC'' package to draw receiver operating characteristic (ROC) curves to evaluate the effectiveness of our model. Finally, the independence of the risk score in HCC was analyzed by Cox regression, and nomogram and calibration curves were plotted. In addition, the model was validated using an ICGC external dataset.

2.10.   Statistical Analysis

In this study, R software (version 4.0.3) was used for data analysis and graphical visualization. The differences in the expression of PAQR4 in HCC were analyzed using the Mann–Whitney U test. Quantitative data were compared and analyzed using the Student's t-test. The Wilcoxon rank sum test or Kruskal–Wallis test, and logistic regression were performed to analyze the associations between PAQR4 and clinicopathological data. Survival differences were analyzed by log-rank test. The independence of variables was analyzed by univariate Cox regression and multivariate Cox regression. A P value < 0.05 was considered as statistically significant.

3.   Results

3.1.   The expression level and prognostic value of PAQR4 in HCC

The TIMER database was used to assess the PAQR4 expression in different cancers. We found that PAQR4 expression was significantly upregulated in lung cancer, gastric cancer, and HCC (Figure 2A). Meanwhile, PAQR4 was highly expressed in HCC tissues compared to adjacent tissues in the TCGA dataset (P < 0.001, Figure 2B, C). Furthermore, we found that PAQR4 was also upregulated in HCC from the ICGC validation dataset (P < 0.001, Figure 2D), GSE14520 validation dataset (P < 0.001, Figure S1A) and GSE36376 validation dataset (P < 0.001, Figure S1B). Moreover, qRT-PCR analysis of six pairs of clinical HCC tissues confirmed that PAQR4 was upregulated in HCC (P < 0.05, Figure 2E). The cell line showed that the expression level of PAQR4 mRNA was higher in hepatoma cells (including, MHCC97H, SKHEP1 and LM3) than in normal hepatocyte LO2 (P < 0.05, Figure 2F). The ROC curve indicated that PAQR4 had a significant diagnostic value for HCC (AUC = 0.907, P < 0.001, Figure 2G). In addition, the prognostic significance of PAQR4 in HCC was analyzed using the Kaplan-Meier survival curves. We concluded that HCC patients with high PAQR4 expression had a poor prognosis in the TCGA dataset (P < 0.05, Figure 3A) and ICGC validation dataset (P < 0.05, Figure 3B).

Figure 2.  The expression level and diagnostic value of PAQR4 in HCC. (A) PAQR4 expression in different tumor tissues and adjacent normal tissues from the TIMER database. (B) The expression level of PAQR4 in HCC tissues (n = 374) and adjacent normal tissues (n = 50) from the TCGA dataset. (C) The expression level of PAQR4 in paired HCC tissues (n = 50) and adjacent normal tissues (n = 50) from the TCGA dataset. (D) The expression level of PAQR4 in HCC tissues (n = 232) and adjacent normal tissues (n = 202) from the ICGC dataset. (E) The relative expression of PAQR4 mRNA in HCC tissues (n = 6) and adjacent normal tissues (n = 6) was analyzed by qRT-PCR. (F) The relative expression of PAQR4 mRNA in hepatoma cells was analyzed by qRT-PCR. (G) The diagnostic value of PAQR4 for HCC in the TCGA dataset. ^*P < 0.05, ^**P < 0.01, ^***P < 0.001, ^****P < 0.0001.

DownLoad: Full-Size Img PowerPoint

Figure 3.  The association between PAQR4 expression and overall survival of patients with HCC. (A) Kaplan-Meier survival analysis of PAQR4 in HCC from the TCGA dataset. (B) Kaplan-Meier survival analysis of PAQR4 in HCC from the ICGC dataset.

DownLoad: Full-Size Img PowerPoint

3.2.   Association between PAQR4 expression and clinical features

A correlation analysis showed that the expression level of PAQR4 was associated with age (P = 0.013), AFP (P < 0.001), histological grade (P < 0.001), T classification (P = 0.005), clinical stage (P = 0.018) and survival status (P = 0.025) in the TCGA dataset (Figure 4). Table 1 shows that PAQR4 expression was also significantly associated with age (P = 0.020), AFP (P = 0.005), histological grade (P < 0.001), T classification (P = 0.006), and clinical stage (P = 0.014) by logistic regression analysis. Furthermore, univariate Cox analysis suggested that clinical stage and PAQR4 were associated with OS of HCC patients in the TCGA dataset (all P < 0.001, Figure 5A). Multivariate Cox regression analysis revealed that clinical stage and PAQR4 were the independent risk factors for HCC (all P < 0.05, Figure 5B).

Figure 4.  Association between PAQR4 expression and clinicopathological parameters. (A) Age. (B) Gender. (C) AFP. (D) Histological grade. (E) T classification. (F) N classification. (G) M classification. (H) Clinical stage. (I) Survival status.

DownLoad: Full-Size Img PowerPoint

Table 1.  Logistic regression of PAQR4 expression and clinical pathological characteristics.

Clinical characteristics	Total (n)	Odds ratio in PAQR4 expression	P value
Age (<65 vs ≥ 65)	372	1.639(1.081−2.496)	0.020
Gender (female vs male)	373	0.758(0.489−1.169)	0.211
AFP (positive vs negative)	279	1.978(1.231−3.199)	0.005
Grade (G3+G4 vs G1+G2)	368	2.672(1.728−4.171)	< 0.001
T classification (T3+T4 vs T1+T2)	370	1.957(1.125−3.190)	0.006
N classification (N1 vs N0)	257	3.072(0.387−62.56)	0.334
M classification (M1 vs M0)	271	1.007(0.119−8.497)	0.994
Clinical stage (Ⅳ+Ⅲ vs Ⅰ+Ⅱ)	349	1.849(1.138−3.035)	0.014

---

 | Show Table

DownLoad: CSV

Figure 5.  The prognostic value of PAQR4 was analyzed in HCC by Cox regression. (A) Univariate Cox regression analysis. (B) multivariate Cox regression analysis.

DownLoad: Full-Size Img PowerPoint

3.3.   Potential mechanism of PAQR4 upregulation in HCC

The potential mechanism of PAQR4 upregulation in HCC was explored from the perspective of genetics. A Pearson correlation analysis showed a significant negative correlation between PAQR4 expression and its DNA methylation in the TCGA dataset (P < 0.001, Figure 6A). Although survival analysis showed no significant difference in the OS of HCC between PAQR4 hypermethylation and PAQR4 hypomethylation (P = 0.227, Figure 6B), PAQR4 hypomethylation & high expression was worse than PAQR4 hypermethylation & low expression (P = 0.018, Figure 6C). The correlation between PAQR4 expression and CNV was explored using the cBioPortal database. We showed that four patients had DNA amplification, 38 patients had DNA copy gain, and PAQR4 copy gain correlated with PAQR4 overexpression (P < 0.05, Figure 6D). Regulator prediction identified hsa-miR-125b-5p as a potential miRNA that regulates PAQR4 expression. Meanwhile, hsa-miR-125b-5p was downregulated in HCC, which was negatively correlated with PAQR4 expression (P < 0.001, Figure 6E, F).

Figure 6.  DNA hypomethylation, DNA copy gain, and miR-125b-5p downregulation contribute to PAQR4 upregulation in HCC. (A) Correlation between PAQR4 methylation and its expression. (B) Survival analysis based on PAQR4 methylation. (C) Survival analysis based on both of the PAQR4 methylation and expression. (D) The expression level of PAQR4 in different CNV groups. (E) The expression level of has-miR-125b-5p in HCC tissues and adjacent normal tissues. (F) Correlation analysis for has-miR-125b-5p expression and PAQR4 expression.

DownLoad: Full-Size Img PowerPoint

3.4.   Analysis of PAQR4 mutation in HCC from cancer genomics

The mutation characteristics of PAQR4 in multiple cancers were analyzed using the cBioPortal database. We found that the frequency of genetic variation in PAQR4 in HCC was less than 1% (Figure 7A). The mutation site of PAQR4 in HCC was shown in Figure 7B. Specifically, the mutant types of PAQR4 were investigated in HCC using the COSMIC database (https://cancer.sanger.ac.uk/cosmic). As shown in Figure 7C, missense mutation was the most frequent type of PAQR4 mutation in HCC (77.78%), followed by synonymous mutations (11.11%). Nucleotide alterations primarily consisted of G > T (87.50%) and C > T (12.50%) mutations (Figure 7D).

Figure 7.  Mutation characteristics of PAQR4 gene in HCC. (A) The mutation frequency of PAQR4 in different cancers from the TCGA dataset. (B) Mutation sites of PAQR4 in HCC from the TCGA dataset. (C) Different mutation types of PAQR4 gene in HCC. (D) Different base mutation proportions of PAQR4 gene in HCC.

DownLoad: Full-Size Img PowerPoint

3.5.   Association of TIICs and immune checkpoints with PAQR4

The correlation between TIICs and PAQR4 was discussed using the TIMER database. The expression level of PAQR4 and the infiltration level of immune cells (including B cells, T cells, dendritic cells, etc.) were significantly positively correlated in HCC (all P < 0.05, Figure 8A). In addition, PAQR4 was positively correlated with immune checkpoint molecules such as PDCD1 (PD1), HAVCR2 (TIM-3), TIGIT, CTLA-4, CD274 (PD-L1) and LAG-3 (all P < 0.05, Figure 8B).

Figure 8.  Association between PAQR4 and TIICs or immune checkpoint molecules. (A)The correlation of TIICs with PAQR4. (B)The correlation of PAQR4 with immune checkpoint molecules.

DownLoad: Full-Size Img PowerPoint

3.6.   Gene set enrichment analysis

The PAQR4-related signaling pathways in HCC were explored using the GSEA. We found that the high PAQR4 expression phenotype was dramatically enriched in the cell cycle, pathways in cancer, homologous recombination, Notch pathway, mTOR pathway, ErbB pathway, Wnt pathway, MAPK pathway, apoptosis and p53 pathway (Figure 9, Table2).

Figure 9.  Enrichment plots from gene set enrichment analysis.

DownLoad: Full-Size Img PowerPoint

Table 2.  The results of gene set enrichment analysis.

Gene set name	NES	NOM p-val	FDR q-val
KEGG_CELL_CYCLE	2.011	0	0.0091
KEGG_NOTCH_SIGNALING_PATHWAY	2.001	0	0.0066
KEGG_HOMOLOGOUS_RECOMBINATION	1.950	0	0.0071
KEGG_MTOR_SIGNALING_PATHWAY	1.919	0	0.0085
KEGG_PATHWAYS_IN_CANCER	1.909	0	0.0070
KEGG_WNT_SIGNALING_PATHWAY	1.904	0	0.0069
KEGG_ERBB_SIGNALING_PATHWAY	1.862	0	0.0082
KEGG_APOPTOSIS	1.855	0	0.0086
KEGG_MAPK_SIGNALING_PATHWAY	1.824	0	0.0109
KEGG_P53_SIGNALING_PATHWAY	1.792	0	0.0136
Note: NES: normalized enrichment score; NOM: nominal; FDR: false discovery rate.

---

 | Show Table

DownLoad: CSV

3.7.   Bioinformatics analyses of the PAQR family

Upon analysis of the PPI network, the results indicated that PAQR family genes have strong protein homology and co-expression (Figure 10A). GO functional annotation analysis revealed that PAQR family genes were markedly enriched in response to hormone, integral component of membrane, molecular transducer activity and other functions (Figure 10B). KEGG pathway enrichment demonstrated that PAQR family genes were enriched in adipocytokine and AMPK pathways (Figure 10B). A Pearson correlation matrix plot of PAQR family genes is shown in Figure10C. The gene expression heatmap of the PAQR family in HCC was shown in Figure 10D.

Figure 10.  Bioinformatics analysis of PAQR family genes. (A) Protein–protein interaction network of 11 PAQR family genes. (B) GO functional annotation and KEGG pathway analysis of the PAQR family. (C) Pearson correlation matrix plot of PAQR family genes. (D) The heatmap visualized of PAQR family genes expression in 374 HCC samples and 50 normal liver samples.

DownLoad: Full-Size Img PowerPoint

3.8.   Establishment and validation of a prognostic mode

To evaluate the prognosis of 11 PAQR family genes in HCC, a univariate Cox regression analysis identified seven genes (PAQR4, PAQR5, PAQR7, PAQR8, PAQR9 and PAQR11) in the TCGA dataset (Figure 11A). Then, multivariate Cox analysis identified three genes (PAQR4, PAQR8 and PAQR9) as the independent risk factors for HCC and included these genes to construct a prognostic model (Figure 11B). The risk formula is constructed through the correlation coefficient (risk score = 0.076614 × PAQR4 expression + 0.070508 × PAQR8 expression + 0.034293 × PAQR9 expression). We utilized the risk formula to calculate the value of the risk score for each patient in the TCGA training set and ICGC test set. All the patients were separated into high-risk and low-risk groups with the median risk value. The risk score, survival status and gene expression cluster plots were drawn in the TCGA training set (Figure 12A−C) and ICGC test set (Figure 12D−F). Kaplan-Meier survival curve analysis demonstrated that patients in the high-risk group had a poor prognosis in the TCGA training set (P < 0.001, Figure 13A) and ICGC test set (P < 0.001, Figure 13B). Meanwhile, the efficacy of the prognostic model was assessed by ROC curve, which showed high accuracy in the TCGA training set (AUC = 0.712, Figure 13C) and ICGC test set (AUC = 0.698, Figure 13D).

Figure 11.  The prognostic value of the PAQR family in HCC by Cox regression analysis. (A) Univariate Cox regression analysis of the PAQR family. (B) The hazard ratios and coefficient of PAQR family prognostic genes were analyzed by multivariable Cox regression.

DownLoad: Full-Size Img PowerPoint

Figure 12.  Risk score, survival status and the gene expression cluster plots of high-risk and low-risk groups. (A) Risk score, (B) survival status, and (C) the gene expression cluster plots of HCC patients in the TCGA training set. (D) Risk score, (E) survival status, and (F) the gene expression cluster plots of HCC patients in the ICGC test set.

DownLoad: Full-Size Img PowerPoint

Figure 13.  Internal and external validation of the prognostic model. (A) Kaplan-Meier survival analysis of high-risk and low-risk groups in the TCGA training set. (B) Kaplan-Meier survival analysis of high-risk and low-risk groups in the ICGC test set. (C) ROC curve to evaluate the performance of the TCGA dataset model. (D) ROC curve to evaluate the performance of the ICGC dataset model.

DownLoad: Full-Size Img PowerPoint

3.9.   Construction of nomogram

The prognostic independence of the risk score was analyzed using the Cox regression in the TCGA dataset. We found that clinical stage and risk score were independent prognostic factors for HCC patients by univariate Cox and multivariate Cox regression analyses (all P < 0.05, Figure 14A, B). The risk score and clinical stage were included to construct a nomogram to predict the survival of patients with HCC (Figure 15A). The calibration curves revealed good agreement between the nomogram-predicted probabilities and actual probabilities in the 1-, 3- and 5-year OS of HCC patients (Figure 15B−D). Meanwhile, areas under the ROC curves of 1-, 3- and 5-year OS were 0.726, 0.709 and 0.687, respectively (Figure 15E). Additionally, the C-index of the prediction model was 0.65, which indicated a satisfactory prediction performance.

Figure 14.  The relationship between the risk score and the prognosis of HCC patients in the TCGA dataset. (A) Univariate Cox analysis of the risk score and clinicopathological parameters. (B) Multivariate Cox analysis of the risk score and clinicopathological parameters.

DownLoad: Full-Size Img PowerPoint

Figure 15.  The nomogram and calibration curve based on the risk score. (A) A nomogram was constructed based on the risk score. (B) The calibration curve of 1-year OS. (C) The calibration curve of 3-year OS. (D) The calibration curve of 5-year OS. (E) ROC curves to assess 1-, 3- and 5-year OS of the nomogram.

DownLoad: Full-Size Img PowerPoint

4.   Discussion

Liver cancer is a highly lethal malignant tumor that seriously affects human health worldwide. An increasing number of studies have identified a variety of key genes affecting the occurrence, progression, and prognosis of HCC ^[28,29]. At present, PAQR4 has been reported to be overexpressed in multiple malignancies and is associated with poor prognosis ^[15,16,30]. In this study, our analysis indicated that PAQR4 was overexpressed in HCC and had a certain diagnostic value for HCC patients. This suggests that PAQR4 may serve as a biomarker for the molecular and pathological diagnosis of HCC. At the same time, patients with high PAQR4 expression had a poor prognosis, which can be used as a potential prognostic indicator for HCC. Previous studies have indicated that PAQR4 is linked to the TNM stage and distant metastasis in lung cancer ^[30]. Notably, this study yielded similar results in that PAQR4 expression was significantly related to age, AFP, grade and clinical stage of HCC patients. Meanwhile, the PAQR4 expression levels increased with increasing clinical stage, suggesting that PAQR4 may be closely related to the occurrence and progression of HCC.

Our study investigated the underlying mechanism of PAQR4 upregulation in HCC from genetic perspective. First, CNV affects the gene transcription or translation, and plays an essential role in tumorigenesis and progression ^[31]. We found that PAQR4 had DNA copy gain and amplification in some patients, which was significantly correlated with PAQR4 upregulation in HCC. Second, DNA methylation is a typical epigenetic modification process. DNA hypomethylation promotes cancer cell proliferation and metastasis, leading to chromosome instability and gene upregulation ^[32,33]. In the present study, we found a significant negative correlation between PAQR4 expression and PAQR4 methylation, suggesting that PAQR4 hypomethylation might be one of the reasons for the upregulation of PAQR4 expression. Finally, miRNAs can promote mRNA degradation and inhibit gene expression, and their aberrant expression can cause tumorigenesis, apoptosis, immune surveillance, and so on ^[34,35]. It has been shown that miR-125b-5p is significantly downregulated in bladder cancer, HCC and breast cancer ^[36,37,38]. We found that miR-125b-5p is downregulated and negatively correlated with PAQR4 expression in HCC, which may be another potential mechanism leading to PAQR4 upregulation. In summary, these results suggest that PAQR4 copy number gain, PAQR4 hypomethylation, and miR-125b-5p downregulation may be possible mechanisms of PAQR4 upregulation in HCC.

Genetic alterations have been related to multiple malignancies ^[39,40]. Additionally, we found that the genetic alteration frequency of PAQR4 was less than 1% in HCC patients in the TCGA dataset. Missense mutations and nucleotide alterations of G > T were the most common types of PAQR4 mutations in HCC. As for tumor immune infiltration cells and immune checkpoints, PAQR4 was significantly correlated with TIICs, PD-1, CTLA-4 and PD-L1. This indicates that PAQR4 is closely related to the immunity of HCC, which might enhance the immune checkpoint inhibitor therapy efforts or even become a new target in HCC immunotherapy.

Moreover, high PAQR4 expression was enriched in some carcinogenesis-related pathways, such as the cell cycle, mTOR, Wnt, MAPK, p53 and other signaling pathways. PAQR4 can reduce the ubiquitination of cyclin-dependent kinase 4, and regulate cell proliferation and tumorigenesis ^[41,42]. Guo et al. reported that the abnormal expression of upstream and downstream genes of mTOR leads to an imbalance in the mTOR pathway, thereby promoting the progression of HCC ^[43]. Meanwhile, both the MAPK and Wnt pathway played crucial roles in multiple processes, and their abnormal sensitization is related to cell atypia and carcinogenesis ^[44,45]. Xu et al. indicated that PAQR4 inhibits the degradation of Nrf2 protein by mediating the Keap1 ubiquitination process, and then promotes chemoresistance in lung cancer ^[17]. Wu et al. determined that PAQR4 enhances lung cancer cell proliferation and metastasis through the CDK4-pRB-E2F1 pathway ^[30]. A recent study reported that PAQR4 promotes prostate cancer cell proliferation, migration, and epithelial mesenchymal transformation (EMT) by mediating the PI3K/Akt pathway ^[15]. Feng et al. showed that PAQR4 promotes cell invasion and EMT in gastric cancer ^[46]. In summary, PAQR4 may promote and influence the carcinogenesis of HCC by regulating various pathways.

In this study, Cox analysis revealed that the three PAQR family genes, including PAQR4, PAQR8, and PAQR9 were independent prognostic factors for HCC. Sinreih et al. reported that PAQR8 is downregulated and is a prognostic biomarker in endometrial carcinoma ^[13]. Del et al. indicated that high PAQR9 expression is associated with poor prognosis in glioblastoma ^[47]. To better evaluate the individualized survival anticipation and plan the short-term follow-up of individual treatment of HCC patients, a prognostic model based on PAQR family genes was constructed using Cox regression analysis. In this study, we included 367 HCC patients from TCGA as the training set and 232 HCC patients from ICGC as the test set. Our model effectively evaluated the patient survival outcomes and was validated in an external dataset, which showed that the model had strong predictive power and value. Furthermore, the nomogram constructed based on the risk score had a good predictive efficiency for 1-, 3- and 5-year survival of HCC patients. We also assessed the accuracy and consistency of the nomogram using the ROC curve and the calibration curve. These results were verified through internal and external validations.

Although our study comprehensively sheds light on the relationship between PAQR4 and HCC and the construction of a remarkably reliable clinical prediction model for HCC, several limitations remain. First, it is not available to obtain a large number of liver tumor samples in the short term, so subsequent work is needed to continue to validate the conclusions. Second, the PAQR4-related signaling pathways are analyzed by bioinformatics, which requires further experimental validation in vitro and in vivo. Third, due to the lack of information related to chemotherapy and targeted therapy in the TCGA database, this is crucial for the prognostic evaluation of patients. Finally, although the prognostic model we construct exhibited excellent predictive value in both the TCGA and ICGC validation cohorts, it failed to be validated in clinical practice.

5.   Conclusions

In conclusion, we explore and find that PAQR4 is overexpressed in HCC, which may be due to PAQR4 hypomethylation, PAQR4 copy gain and miR-125b-5p downregulation. HCC patients with high PAQR4 expression have a poor prognosis. Moreover, PAQR4 is closely associated with tumor-infiltrating immune cells, immune checkpoint molecules, and some oncogenic-related pathways in HCC. Our nomogram model could accurately assess the clinical prognostic significance in patients with HCC. More importantly, our study comprehensively and deeply explores the relationship between PAQR4 and HCC and constructs a prognosis model based on the PAQR family, which will provide more scholars with new ideas and methods. In the future, we believe that PAQR4 may serve as a target for HCC immunotherapy and work together to give our prognostic model to clinical workers a better assessment and planning of patient individualized survival anticipation and short-term follow-up.

Acknowledgments

This study was supported by the National Natural Science Foundation of China (No. 31770537); Gansu international science and technology cooperation project (No.18YF1WA113); Major science and technology projects in Gansu Province (No. 20ZD7FA003).

Conflicts of interest

The authors declare that they have no competing interests.

Grant Disclosures

The following grant information was disclosed by the authors: This study was approved by the Institutional Research Ethics Committees of the Lanzhou University Second Hospital (protocol number: 2019A-321).

Appendix

Figure S1.  The expression level of PAQR4 in HCC. (A) The expression level of PAQR4 in HCC tissues (n = 225) and adjacent normal tissues (n = 220) from the GSE14520 dataset. (B) The expression level of PAQR4 in HCC tissues (n = 240) and adjacent normal tissues (n = 193) from the GSE36376 dataset.

DownLoad: Full-Size Img PowerPoint