PCK2:A potential Therapeutic Target and Biomarker for Hepatocellular Carcinoma Based on Weighted Gene Co-Expression Network Analysis

Background Hepatocellular carcinoma(HCC) is the sixth most common cancer worldwide, and is associated with a high mortality rate and poor treatment efficacy. Therefore, it is important to explore the underlying mechanisms and screen for novel biomarkers and therapeutic targets. HCC and non-tumor microarray data was obtained from the Gene Expression Omnibus(GEO) database and 1,612 differentially expressed genes (DEGs) were identified. The identified DEGs were then further examined using a weighted gene co-expression network analysis (WGCNA). Modules were identified and correlations with clinical parameters were examined to ultimately identify a single highly correlated module containing a single hub gene, phosphoenolpyruvate carboxykinase 2 (PCK2).The genes co-expressed with PCK2 were then examined via Gene Ontology(GO) analysis and Kyoto Encyclopedia of Genes and Genomes(KEGG)pathway analysis and were found to be primarily enriched in the area of energy metabolism. The association between PCK2 expression and prognosis and tumor stage was then further verified by using a HCC dataset obtained from The Cancer Genome Atlas (TCGA) database. Ex vivo examination of HCC and adjacent non-tumor samples via quantitative real-time PCR demonstrated that PCK2 is downregulated in HCC. Additionally, HCC cell lines overexpressing PCK2 exhibited a time-dependent reduction in cellular viability. These findings demonstrate that PCK2 is downregulated in HCC and is associated with HCC prognosis and tumorstage, and that PCK2 may affect HCC progression via energy metabolism. this study may help to further elucidate possible HCC progression mechanisms and the findings presented here suggest that PCK2 may serve as a potential HCC therapeutic target or biomarker. The significance of PCK2 in HCC was further verified by examining clinical quantitative real-time PCR (qPCR) data obtained from The Cancer Genome Atlas (TCGA) database. These findings were then further confirmed by evaluating a TCGA dataset. PCK2 and its co-expressed genes were then functionally examined by performing GO and KEGG enrichment analyses, which showed a strong enrichment in metabolic related biological processes. In tumors, metabolic reprogramming is an important feature, with the

3 World Health Organization estimates that by 2030, the number of HCC-related deaths will exceed one million [3]. Additionally, HCC is the second-most lethal tumor after pancreatic cancer and is associated with a five-year survival of 18% [4]. Despitethe availability of many HCC treatments, including a local resection, transplantation, chemotherapy, orimmunotherapy, the overall prognosis remains poor [5].To improve treatment, an understanding of HCC occurrence and development at the molecular level is necessary to identify novel treatments and biomarkers.
Various platforms, including microarray and sequencing based approaches, can be utilized to elucidate the molecular mechanisms associated with tumor occurrence and development and aid in the isolation of potential therapeutic targets and biomarkers. Gene co-expression networks (GCNs) are constructed by utilizing high-throughput gene expression datasets, and aid in the identification of networks associated with certain biological processes for a given disease [6]. Weighted gene coexpression net work analysis(WGCNA) is a powerful analytical tool that explores modular information based on weighted correlation net works and reveals correlations between modules and clinical parameters [7]. By using WGCNA, co-expressed genes can be systematically clustered in a more reliable way and can be combined with corresponding clinical parameters to isolate specific biological processes that are altered in patients with HCC, as well as to discover candidate biomarkers and/or therapeutic targets [8]. WGCNA has been used to study many complex diseases, including breast cancer, prostate cancer, osteosarcoma, and coronary heart disease [9][10][11][12]. While a large number of genes associated with HCC occurrence and development have been identified [13,14], more HCCrelated genes still need to be discovered.
In the present study, differential expression analysis was performed using microarray data obtained from the Gene Expression Omnibus (GEO) database. Differentially expressed genes(DEGs) were analyzed via WGCNA to identify gene modules with HCC biological significance and to discover hub genes for a given module. The results identified a single module with a single hub gene, phosphoenolpyruvate carboxykinase 2(PCK2), which was found to be associated with HCC survival and tumor stage. The significance of PCK2 in HCC was further verified by examining clinical quantitative real-time PCR (qPCR) data obtained from The Cancer Genome Atlas (TCGA) database. 4 Additionally, the effect of PCK2 was examined in vitro by constructing a cell line overexpressing PCK2.
The purpose of the present study was to identify genes associated with HCC occurrence, development or prognosis that may serve as potential biomarkers or therapeutic targets to aid in HCC detection and treatment.

Data collection, processing and differential expression analysis
Gene microarray expression data(GSE76427), and associated clinical data, were obtained from the National Center of Biotechnology Information(NCBI) GEO database(http://www.ncbi.nlm.nih.gov/geo).
The GSE76427 dataset contains complete expressional data for 115 HCC samples and 52 non-tumor samples. The data was downloaded and analyzed using R software v. 3.6.1 www.r-project.org/ .

Weighted gene co-expression network construction
DEGs were utilized to construct the co-expression networks using the WGCNA package in R. Pairwise correlation analyses were performed and an adjacency matrix was calculated using an appropriates soft threshold to generate a scale-free network using the integrated function (pick Soft Threshold) in the WGCNA package. For the clustering analysis, aminimum module size=30(genegroup) was utilized for the gene dendrogramand, and a height cut off=0.25 was used for the module dendrogram.

Identification of clinically relevant modules and candidate hub genes
To identify clinically relevant modules, correlations between the clinical data, including age, sex, tumor-node-metastasis(TNM)staging classification, overall survival (OS) status, and OS time, and then the identified modules was examined.Two methods were used to test the correlation between each module and the development of HCC. Modules with a higher gene significance (GS) or higher correlation of module eigengene(ME) and GS were considered to be more related to the development of the disease and were selected for further analysis. Next, Hub genes were identified based on module connectivity (gene Module Membership>0.8) and clinical trait relationship(gene Trait Significance>0.2).

Correlation analysis between hub genes and survival and tumor TNM stage
The optimal cut-off value was determined by the surv_cutpoint function in survminer (R package).
Based on the cutoff value, patients were divided into high expression or low expression groups and a Kaplan-Meier survival curve was plotted. To determine possible correlations between hub genes and TNM stage, a linear regression analysis was performed.

Functional enrichment analysis
To further characterize the underlying mechanisms of hub gene, the co-expression genes of hubgene were screened and Gene Ontology(GO) annotation and Kyoto Encyclopedia of Genes and Genomes(KEGG) pathway analysis were performed using Database for Annotation, Visualization and Integrated Discovery(DAVID; http://david.abcc.ncifcrf.gov/).

TCGA validation dataset
To further validate the findings based on the GEO dataset, a HCC dataset was downloaded from the TCGA database, with the associated clinical data obtained from UCSC XENA(https://xena.ucsc.edu/). A total of 343 samples were included in the analysis, with some samples removed due to incomplete survival and/or TNM stage information. The 343 samples were then divided into high expression and low expression groups based on the median expression of the hub gene. Then, a Kaplan-Meier survival curve and a regression analysis were performed to determine correlations between gene expression and survival or TNM stage.

Validation of the hub gene via quantitative real-time PCR (qPCR)
A total of 24 HCC patients were recruited to obtaining tumor and adjacent non-tumor tissues. The patients provided informed written consent and this study was approved by the Ethics Committee of

Construction of cell lines over-expressing PCK2 and cell-proliferation assay
SMMC7721 and Huh-7 cell lines were obtained from the College of Basic Medical Sciences, JilinUniversity, Changchun, Jilin, China.The cells were cultured in DMEM medium supplemented with 10% fetal bovine serumand were incubated at 37ºC with 5% CO 2 . The cells were collected at the logarithmic phase and were seeded into a 96-well plate at 5,000 cells/well.Each cell line was then transfected with PCK2 overexpression plasmids(OriGeneTechnologies,Inc.China) or with empty plasmids using Trans Intro EL Transfection Reagent(Transgen,China) according to the manufacturer's protocols.Successful transfection was verified by examining PCK2 expression via qPCR as described above.
To evaluate the effect of PCK2 overexpression on cellular proliferation, a MTT assay was performed. At 24, 48, and 72h post-transfection, 100μl of MTT solution (0.5mg/ml)was added to each well. The samples were then incubated for 4h at 37°C, with the solution then discarded and replaced with 150μl of dimethylsulfoxide (DMSO).The absorbance was measured at 490nm using a plate 7 reader.

Statistical analysis
Data obtained from the qPCR and MTT analyses were analyzed using GraphPad Prism (version7.0),with differences among groups statistically evaluated by using a two-tailed Student's ttests. Differences were considered statistically significant at P<0.05, with all samples examined in triplicate.

Weighted gene-coexpression network construction
After removing any outliers, the 115 HCC samples, which comprised 1,612 DEGs, were clustered and possible correlations with clinical data were evaluated (Figure2(A)).Based on the criterion for a scalefree topology, a soft threshold power of 4 was chosen,with a scale-free topology index(R 2 )=0.85 obtained (Figure2(B)).

Clinical trait associated modules and hub gene screening
After applying the soft threshold, eight modules were screend out for further analysis. Of the screened modules, the brown module( Figure 3)had the highest GS thus the highest correlation with HCC development than did the other gene modules (Figure3 (A-C)).Therefore,the brown module was selected for further analysis.When evaluating the module connectivity(gene Module Membership>0. 8) and clinical trait relationship(geneTraitSignificance>0.2),only PCK2 was isolated as a hub gene (Figure3(D)), it was downregulated in HCC tissues in the dataset. 8

Survival and tumor stage correlation analysis of the hub gene
For survival analysis in association with the hub gene PCK2, an optimal cutoff value was identified according to survminer, and the patient data were divided into high expression and the low expression groups. Kaplan-Meier survival curves were plotted, and the result showed that the survival rate of the high expression group was significantly reduced when compared to the low expression group(P<0.05 ;Figure4(A)).To determine a potential association between PCK2 and different TNM stages, a linear regression analysis was performed, with P< 0.05 indicating a significant correlation (Figure4(B)).

Functional enrichment analysis of genes co-expressed with PCK2
Genes co-expressed with PCK2 were screened with a significance threshold value of P<0.01 and a correlation coefficient>0.45. A total of 250 genes were utilized for enrichment analysis. GO enrichment analysis was performed, and the top three biological processes included the smallmolecule catabolic process, organic acidcatabolic process,and the carboxylic acid catabolic process ( Figure 5(A)).KEGG pathway analysis was also preformed, with the top three enriched pathways relating to peroxisomes, fatty-acid degradation, and valine,leucine and isoleucine degradation ( Figure   5(B)).

Validation of PCK2 using a TCGA dataset
The relationship between PCK2 expression level and HCC prognosis was further verified by utilizing gene expression profiles and associated clinical data for 343 HCC patients from the TCGA database.The results shown PCK2 was downregulated in HCC and low expression has a worse prognosis, which were consistent with the previous results(Figure6).

Ex vivo and in vitro validations via qPCR and MTT assay
To further explore the role of PCK2, ex vivo verification was performed by examining 24 HCC tumor and adjacent non-tumor patient samples. In the HCC tumor samples, PCK2 was significantly downregulated relative to the matched non-tumor samples(P<0.05; Figure7(A)).
The impact of PCK2 expression was also examined in vitro by constructing SMMC7721 and Huh-7 cell overexpressing PCK2. Successful transfection was confirmed at 24h post-transfection using qPCR (Figure7(B)). To assess cell viability in the presence of PCK2 overexpression, MTT assays were performed at 24, 48, and 72h post-transfection. The results showed that PCK2 overexpression reduces cell viability in a time-dependent manner after 48h post-transfection (Figure7(C)).

Discussion
HCC is a devastating disease worldwide. Although multivariate analyses of HCC development have been extensively documented, the molecular mechanisms underlying HCC remain poorly understood [15]. Therefore, elucidating disease-specific gene expression patterns has helped to uncover some of the underlying mechanisms of different diseases [16].When compared to other analytical approaches, WGCNA provides a higher reliability and biological significance due to both genetic modules and clinical parameters being simultaneously examined [17]. Furthermore, WGCNA has been widely used to screen genes associated with pathogenesis or cancer clinical parameters, including lung cancer [18],colon cancer [19],melanoma [20],and osteosarcoma [21]. In the present study, the mainaim was to utilize WGCNA to aid in the identification of key genes associated with HCC development and prognosis. Any findings would then be preliminarily verified at the tissue and cellular level to ultimately aid in the identification of potential HCC biomarkers or therapeutic targets.
When analyzing microarray data obtained from a GEO dataset (GSE76427), 1,612 HCC DEGs were identified when compared to the obtained non-tumor samples. Co-expressional analysis (WGCNA) was then performed and genes were consolidated into modules. Of the identified modules, 8 modules were determined to be associated with the examined clinical parameters, including age, sex, tumor stage, OS status and OS duration. The brown module was found to be the most clinically significant, with PCK2 confirmed as the module hub gene. These findings were then further confirmed by evaluating a TCGA dataset. PCK2 and its co-expressed genes were then functionally examined by performing GO and KEGG enrichment analyses, which showed a strong enrichment in metabolic related biological processes. In tumors, metabolic reprogramming is an important feature, with the metabolism playing an important role in tumor anti-oxidative stress, energy production and biosynthesis [22].Therefore, key metabolic genes may play an important role in regulating metabolic networks and may affect various biological processes associated with tumors. Furthermore, PCK2 overexpression in HCC cell lines was shown to reduce cell proliferation, thus further suggesting that PCK2 may be a key gene that affects HCC occurrence and development.
There are two PCK genes, cytoplasmic PCK1 and mitochondrial PCK2. In the cells that repopulate melanoma tumor cells, PCK2 has been shown to act as hub between tricarboxylic acid (TCA) cycle, gluconeogenesis and glycolysisby modulating the conversion of mitochondrial oxaloacetate (OAA) to phosphoenolpyruvate and by regulating the glucose carbon flow [19]. Furthermore, PCK2 downregulation was shown to promote biosynthesisand the transportation of citricacid from the mitochondria to the cytosol; while overexpression hindered cell growth and tumorigenesis [19], which is consistent with the findings presented herein.
Phosphoenolpyruvate carboxy kinase (PEPCK), which is encoded by PCK, has two isoenzymes, cytosolic PEPCK1(PEPCK-C) and mitochondrial PEPCK2(PEPCK-M), which are encoded by PCK1 and PCK2, respectively. While both isoforms are expressed in the liver, PEPCK-M is also expressed in tumors and a variety of cell types [20]. In the mitochondria, PEPCK-M catalyzesthe conversion of oxaloacetatetophosphoenolpyruvate(PEP) and plays a role in gluconeogenesis, while also being present in a variety of non-gluconeogeni c tissues. Furthermore, in MCF-7 mammary carcinomacells, PEPCK-M knockdown in the presence of stress, such as amino-aciddeprivation or ER stress, reduced cell growth [23]. In another study examining non-small-cell lung carcinomas and adjacent non-tumor tissues, PEPCK activity was three times higher in the tumor samples when compared to the nontumor. Moreover, when inhibiting PCK2 via 3-mercaptopicolinate or small interfering RNAs in various cell lines, glucose depletion-induced apoptosis was significantly enhanced in A549 and H23 cells, but not in H1299 cells. Hence, PCK2 is actively expressed in lung cancer and may aid in cancer cells adapting to low-glucose conditions [24]. Moreover, Zhao et al. found that in prostate cancer, PCK2 is associated with proliferation and metastasis, with higher PCK2 expression levels associated with more aggressive tumors and a lower survival rate [25].Overall, these findings are consistent with the findings presented herein that also show that a high PCK2 expression levels are associated with a lower survival rate.
Little has been reported regarding the expression pattern and function of PCK2 in HCC. In the present study, PCK2 expression was downregulated in HCC tissues and played an inhibitory role in HCC cell lines, but the underlying mechanisms were not elucidated. As a key gluconeogenesis enzyme, the main function of PEPCK is cataplerosis [26], they may play synergistic role with various proteins.
Enrichment analysis suggested that co-expression genes of PCK2 were enriched in energy metabolism. In addition, as a similar structure and different subcellular localization with PEPCK-M, PEPCK-C has been reported to severe in cataplerosis in HCC cells under glucose deficient conditions, with PEPCK-C overexpression significantly reducing glycolysis and the TCA cycle, while increasing reactive oxygen species(ROS) levels and subsequently leading to oxidative stress. Furthermore, PEPCK-C overexpression induces apoptosis and inhibits liver tumorigenesis in vivo [27,28], it is similar to our study on PCK2 overexpression in HCC cells. Thus, it is possible that PEPCK-M overexpression may hinder glycolysis and the TCA cycle too, thus leading to an intracellular energy imbalance that could inhibit cell survival. However, further experimentation is required to confirm this hypothesis.

Conclusion
In the present study, microarray data and clinical information were analyzed and integrated via WGCNA, and PCK2 was identified as a key gene that is highly correlated with HCC survival and tumor stage,with these findings further confirmed by utilizing a TCGA dataset.Furthermore, when examining HCC cell lines overexpressing PCK2, PCK2 expression weakened in a time-dependent manner.While the underlying mechanism for this remains unclear, it might be due to PCK2-mediated prevention of glycolysisand the TCA cycle. However, further experiments are needed to verify or refute this hypothesis.Taken together, PCK2 maybe able to serve as a potential HCC biomarker and/or therapeutic target, and may provide a new strategy for HCC diagnosis and therapy.

Ethics approval and consent to participate
This study was approved by the Research Ethics Committee of Sencond Hosipital of Jilin University.

Consent for publication
All participants gave written consent for the publication of the identifying photographs.