Glia Maturation Factor Beta as a Novel Biomarker and Therapeutic Target for Hepatocellular Carcinoma

Hepatocellular carcinoma (HCC) is one of the most common types of cancer. The novel sensitive biomarkers and therapeutic targets are urgently needed for the early diagnosis of HCC and improvement of clinical outcomes. Glia maturation factor-β (GMFB) is a growth and differentiation factor for both glia and neurons and has been found to be tightly involved in inflammation and neurodegeneration conditions. In our study, the expression level of GMFB was significantly up-regulated in patients with HCC and positively co-expression with tumor node metastases (TNM) stage and histopathological grade of HCC. The high expression level of GMFB was remarkably associated with poor overall survival, which mainly occurred in males rather than females. Multivariate analysis revealed GMFB to be an independent prognostic factor for overall survival in patients with HCC. Results of Gene Ontology (GO) and KEGG pathways analysis showed that down-regulation of pathways related to protein translation and mitochondria function were enriched. Protein-protein interaction analysis revealed the central role of mitochondria protein in HCC. The downregulation of genes involved in glycolysis and gluconeogenesis was observed among the co-expression genes of GMFB. Knockdown of GMFB in Hep3B significantly inhibited proliferation, migration, and invasion of Hep3B cells, and also downregulated the expression levels of some of metal matrix proteinase (MMP), increased mtDNA copy number and loss of mitochondrial transmembrane potential. GMFB influences the malignancy rate of HCC possibly through regulation of the expression of MMPs, mtDNA function and glycolysis. We proposed that GMFB was a promising HCC diagnostic and prognostic biomarker and therapeutic target in HCC.


INTRODUCTION
Hepatocellular carcinoma (HCC) is the main component of primary liver cancer, accounting for 75%-85%, and remains the second leading cause of the death rate of cancers worldwide (1,2). The major risk factors of HCC showed regional differences. Among the most important risk factors, chronic infection with hepatitis B virus (HBV) or hepatitis C virus (HCV), heavy alcohol intake, obesity and type 2 diabetes were highlighted (3). Although HCC can be prevented by the vaccine against HBV, the incidence and mortality of HCC have been increasing rapidly for the past several years worldwide and represents a considerable public health burden (4)(5)(6)(7)(8).
Molecularly targeted therapeutics are important methods for treating advanced HCC (15). To date, a large number of biomarkers have been served as therapeutic targets for HCC. For instance, as a multi-targeted tyrosine kinase receptor inhibitor Sorafenib exerts its function by targeting CRAF、 BRAF、VEGFRl-3、PDGFR-p、cKIT、FLT-3 and RET (16). Although remarkable progress has been achieved, there are no robust biomarkers that could be applied in the early diagnosis and treatment of HCC (17). However, a limited effect with a fiveyear survival rate of 18% in HCC patients was obtained (18,19). Therefore, there is an urgent need to develop a novel, highly sensitive and specific biomarker for the diagnosis, prognosis and therapeutic target for HCC (20). There is also a remarkable gender disparity of HCC incidence, primarily males (21). The current study provides evidence that androgen in males could play a critical role in hepatocarcinogenesis and development, and estrogen could protect women from hepatocarcinogenesis (22,23). Therefore, the gender-related differences should be given special consideration in individual HCC treatment from the perspective of precision medicine.
GMFB is a growth and differentiation factor expressed predominantly in the central nervous system (CNS) and testis (24,25). As for liver, the expression of GMFB was confirmed by previous proteomics analysis (26,27). GMFB is involved in growth and differentiation in the vertebrate brain, neutrophil chemotaxis, migration of monocyte, migration, and adherence of T lymphocytes (28)(29)(30)(31). GMFB was found to be upregulated in several neuroinflammation and neurodegeneration conditions (32). Increased GMFB expression was found highly associated with multiple types of cancer, including glioma and ovarian cancer, and GMFB overexpression was reported to be coexpression with poor prognosis in ovarian cancer and glioma (33,34). To date, there have been no studies of the effects of GMFB on HCC, therefore its roles in HCC progression remain unclear.
This study aimed to characterize the association between GMFB and HCC, and reveal a potential underlying mechanism of GMFB of HCC. To verify analysis results, a series of experiments were carried out. The expression level of GMFB affects mtDNA copy number, mitochondrial membrane potential and MMPs expression levels in Hep3B cell, and impairs cell migration, invasion and adhesion. To our knowledge, this is the first study to report the diagnostic and prognostic values of GMFB in HCC and analyze the genderbased correlation between GMFB and overall survival in patients with HCC. Our study provided a new potential diagnostic and prognostic marker, and a novel bio-target for the treatment of HCC.

Ethics Statement
Our study was approved by the Academic Committee of Tongji University and conducted in accordance with the Declaration of Helsinki. We obtained all the data from the online databases, therefore the informed consent for data collection had already been obtained.

Oncomine Database Analysis
ONCOMINE (http://www.oncomine.org) is an integrated translational bioinformatics platform composed of datasets and analyses (35). Datasets include samples represented as microarray data measuring either mRNA expression or DNA copy number, which is used to set up analyses on groups of interest like cancer versus normal. Oncomine analyses result from computations that are performed on samples within a dataset. To investigate the transcription levels of GMFB in different types of cancers the Oncomine database was used. Student´s t-tests were performed, with results filtered by the cut-off of p-value <0.05, fold change of >1.5 and gene rank in the top 10%.

UALCAN Database Analysis
UALCAN (http://ualcan.path.uab.edu) is a comprehensive web resource based on 3 RNA-seq databases and the clinical data of 31 cancer types from the The Cancer Genome Atlas (TCGA) database (project ID, TCGA-LIHC) (36). To analyze the levels of the ADF family protein in liver hepatocellular carcinoma. We analyze the transcriptional expression patterns in normal liver tissue and liver hepatocellular carcinoma samples through UALCAN.

Human Protein Atlas Analysis
HPA (https://www.Proteinatlas.org/) is a valuable tool provided immunostaining on tissues and cell lines as well as differential expression analysis of proteins in normal and tumor tissues (37). In this study, we checked the protein expression of GMFB in the HPA database and analyzed the immunohistochemical results of GMFB in tumor tissues and normal tissues (Antibody: HPA053669). cBioPortal cBioportal (http://www.cbioportal.org/) is an open access powerful tool offered to explore, visualize and analyze multidimensional cancer genomics data. The genomic data include somatic mutations, DNA copy-number alterations (CNAs), mRNA and microRNA (miRNA) expression, DNA methylation, protein abundance and phosphoprotein abundance (38). The Oncoprint module in cBioportal was used to provide a graphic summary of major genetic alterations and changes in gene expression. Increases or decreases in mRNA level were based on a Z-score threshold of 2.0 or more standard deviations from the mean of the reference population. The reference population was the samples that are scored diploid for each gene. In our study, the relationship between GMFB in HCC patients and alteration frequency were analyzed through cBioportal. Venn diagram was created using online Venn software (http://bioinformatics.psb.ugent.be/webtools/Venn/).

LinkedOmics and GSE
The LinkedOmics database (http://www.linkedomics.org) is an open-access online biometrics platform that contains 32 cancer types comes from 11,158 patients from TCGA database (project ID, TCGA-LIHC) (39). In this study, we determined the GMFB associated gene enrichment using the "LinkInterpreter" module, which performs enrichment analysis based on Gene Ontology, biological pathways, network modules, among other functional categories.

Metascape, DAVID and Network Analyst Analysis
Metascape and DAVID are online software for gene annotation and gene set enrichment analysis. In this study, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of GMFB co-expression genes were performed using Metascape and DAVID. And the PPI networks were conducted by Metascape. In Metascape, the min overlap was set as 3, the P-value cutoff was 0.01, and the minimum enrichment factor was 1.5. The cutoff value was set as FDR of<0.05. NetworkAnalyst (http://www.networkanalyst. ca) was used to identify the differential expression of GMFB between HCC and normal samples from TCGA database (project ID, TCGA-LIHC) (40).

Construction of Prognostic Models and Survival Analysis
The Kaplan-Meier Plotter database (http://kmplot.com/analysis/) is an online public database, to draw survival plots using publicly available data (41). In this work, our survival analysis was carried on the Kaplan-Meier Plotter database contains survival information for 364 patients with HCC, 250 male patients and 121 female patients with the "Auto select best cutoff" option. The expression levels of GMFB and clinical features were merged to find independent prognostic factors through univariate and multivariate Cox regression analysis. The receiver operating characteristics curve (ROC) analysis was performed using the 'survival ROC' package in R. The area under the ROC curve (AUC) was measured for the prediction of GMFB.

Cell Lines
Hep3B cells purchased from ATCC were used in this study. All cells were maintained in High Glucose Dulbecco's Modified Eagle Medium (DMEM) supplemented with 10% FBS (Gibco, Thermo Fisher Scientific, Spain) at 37°C and 5% CO 2 . and cultured.

Cell Transfection
Hep3B cells were transfected using Lipofectamine 2000 (Invitrogen, USA), according to the manufacturer's instructions. SiRNA and Lipofectamine 2000 were diluted in FBS free DMEM(Gibco) medium and incubated for 5 min, separately. Then mixed together, and incubated for 10-15 min at room temperature (RT) to form the DNA-Lipofectamine complexes. Finally added to the culture plates.

Western Blot Analysis
The samples were lysed with cell lysis buffer (RIPA Lysis Buffer; CAT: P0013B; Beyotime Institute of Biotechnology) to extract whole-cell protein. Protein concentration was quantified using a BCA kit (CAT: P0012; Beyotime Institute of Biotechnology Jiangsu, China), and 20-50 mg of each protein were separated by SDS-PAGE using a 12.5% SDS-PAGE gel. After transfer to a polyvinylidene fluoride membrane (Millipore, USA), membranes were blocked with 5% skim milk at room temperature for 1 hr and incubated overnight at 4°C with primary antibody against GMFB (1:2000; Cat no:10690-1-AP; Proteintech) and ACTB (1:5000; Cat no: 20536-1-AP; Proteintech). Secondary antibodies used were rabbit anti-goat Ig-HRP (1:5000; Proteintech). The membranes were then incubated with secondary antibodies for 2 hrs. Proteins were visualized using ECL (Millipore). The density of the bands was determined using ImageJ software (USA).

Wound Healing Assay
Hep3B cells were seeded in 6-well plates and grown to 90% confluence. The confluent cell monolayers were scratched using a sterile 200µl pipette tip and rinsed with PBS to remove scratched cells. Wound closure was observed for 36 hours. The wound closure rate was calculated as follows: wound closure (%) = (area of initial wound-area of final wound)/area of initial wound ×100.

Transwell Migration and Invasion Assays
For cell invasion analysis, transwell chambers were coated with Matrigel, and the cell migration (without matrigel) assay was also performed. Cell suspensions in FBS free medium were added to upper transwell chambers (24-well, 8-mm pores; BD Labware, USA), While 600 µL DMEM medium supplemented with 10% FBS was added into lower transwell chambers. Cell counting was performed using the cell counter plugin of ImageJ. The migratory capacity of Hep3B cells under various treatments was evaluated through scratch assay.

Measurement of mtDNA Copy Number
Total DNA (mtDNA and gDNA) was isolated using TIANamp Genomic DNA Kit from TIANGEN (TIANGEN, Beijing, China) according to the manual. Relative mtDNA copy number was quantified by quantitative PCR (qPCR) (42). The quantitative PCR was performed by using SuperReal PreMix Plus (SYBR Green) (TIANGEN, Beijing, China) on a BIO-RAD CFX96 Touch Real-Time PCR Detection System. The mtDNA (Dloop, MT-TL1 and MT-ND1) was normalized to nuclear DNA (NCOA3). The primers and primer sequences of mtDNAspecific primers used were as previously reported (43). The relative fold changes were calculated by 2 −DDCt .

Quantitative Real-Time PCR
Total RNA was isolated using the TRIzol reagent (Takara, Dalian, China) strictly according to the instructions. Then RNA was reverse transcribed to cDNA using PrimeScript RT polymerase (Takara, Dalian, China). SYBR Green Master Mix (Tiangen Biotech, China) on a LightCycler 96 Detection System (Roche) was used for RT-qPCR. Data were analyzed using the 2 -DDCt method. Oligos were synthesized by Sangon. Primer sequences were obtained from Origene.

Mitochondrial Membrane Potential
Hep3B cells were seeded onto coverslip placed in 24-well plates and cultured at 37°C in a 5% CO2 incubator. After 24 h, cells were transfected with scramble control siRNA or GMFB siRNA with Lipo2000 (Invitrogen). Scramble siRNA was used as a control. JC1 staining was performed 48 h after transfection. Cells were washed three times with PBS, and then stained with the 2 mM JC-1 dye (Beyotime, Shanghai, China) in the dark for 20 minutes. Cell visualization was performed by a fluorescence microscope.

Statistical Analysis
Image analysis was performed using Image J. Student t-test and one-way ANOVA were used in the statistical analysis. Data were expressed as means ± SEM. The p-value less than 0.05 was considered statistically significant.

Increased Expression of GMFB in Liver HCC
We initially used the Oncomine database to analyze the expression profiles, and found that the GMFB transcription levels were markedly higher in breast cancer, head and neck cancer and liver cancer compared to the normal tissues ( Figure 1A). To study the changes of GMFB mRNA expression levels between HCC and normal liver tissues, we then used the Oncomine and UALCAN database to analyze the differential. The mining of publicly available databases from Roessler Liver 2 Statistics showed that GMFB was up-regulated in HCC ( Figure 1B). Roessler Liver 2 Statistics was containing 12,624 measured mRNAs from 445 samples. Samples were hybridized to Affymetrix Human Genome HT U133A Array. GEO sample IDs were given in Supplementary Table 1. Then, we explored the expression of GMFB in HCC using the UALCAN based on the data resources of The Cancer Genome Atlas database. The results were nearly identical ( Figure 1C). Furthermore, the analysis of TCGA datasets revealed that the expression of GMFB in HCC was associated with TNM stage 1-3 ( Figure 1D) and pathological grade 1-3 ( Figure 1E). Considering the small number of samples of stage 4 and grade 4, the statistical conclusion may not be accurate, and future larger cohort validation is needed. The analysis of TCGA datasets revealed that GMFB was highly expressed in both male and female HCC patients using the UALCAN. And there was no difference was found between male and female HCC patients ( Figure 1F). To solve the imbalance among normal, male and female data in Figure 1F, we downloaded normal (liver), male (HCC) and female (HCC) gene expression data from TCGA further validation. Each sample was assigned a number and 30 samples of each group were randomly selected (Supplementary Table 2). The randomization was computer-based, generated in R (version R 3.4.4). The result indicated that there was a significantly higher expression of GMFB in HCC patients in both genders, and no significant gender differences were found in normal (liver) or HCC tissues ( Figure 1G). Furthermore, Detection of GMFB protein expression in normal and HCC liver tissue of females and males by immunohistochemistry retrieved from the Human Protein Atlas also indicated higher protein expression patterns of GMFB in HCC tissues when compared to normal samples ( Figure 2A). Furthermore, we also performed GMFB copy number variation (CNV) and mRNA levels among normal and HCC tumor samples using Oncoprint algorithm from cBioPortal, and the result showed that 5% of the HCC cases had undergone genetic changes ( Figure 2B).

The Significant Correlation Between Increased GMFB Expression and Poor Overall Survival in HCC
To clarify the associations between the expression levels of GMFB and patients' clinical outcomes, Kaplan-Meier (KM) Plotter survival analysis was performed. The survival curves indicated that high GMFB expression was significantly associated with poor overall survival (OS) of HCC patients (P=0.0042) ( Figure 3A). Next, we analyzed the effects of gender differences on the association between GMFB and OS. A significant difference was observed in male patients ( Figure 3B), but not in females (P=0.36) ( Figure 3C).

GMFB Being an Independent Risk Factor for HCC
The AUC values with 95% CI of GMFB was 0.852, CI in a range from 0.809 to 0.895, indicating a significantly diagnostic accuracy for HCC ( Figure 3D). Univariate analysis by Cox proportional hazard model showed that GMFB, pathologic stage and tumor status were responsible for OS in HCC. Furthermore, the multivariate Cox regression analysis confirmed that GMFB (P=0.039), pathologic stage(P<0.001) and tumor status(P=0.003) were independent factors for predicting the prognosis of HCC patients ( Figure 3E). The clinical characteristics of the HCC patients were shown in Supplementary Table 3.

The Analysis of Co-Expression Genes With GMFB in HCC Patients
To deep mine the underlying mechanisms of the role of GMFB in HCC, LinkedOmics was utilized to analyze mRNA sequencing information from HCC patients in the TCGA. The plot showed the distribution of significant positive correlation with GMFB (red dots) or significant negative correlation (blue dots) in HCC patients ( Figure 4A Table 4). GO (BP, MF, and CC) and KEGG enrichment analyses of co-expression genes were performed. As shown in Figure 4D, we found that the down-regulated co-expression genes of GMFB were associated with the mitochondrial electron transport chain and translation initiation, and the up-regulated co-expression genes of GMFB were relevant to transcription, DNA-template (BP), nucleoplasm (CC), protein binding (MF) and pathway in cancer (KEGG). The detailed information for the co-expression genes was obtained in Supplementary Tables 5, 6. Furthermore, results of Kaplan Meier survival analysis (males and females were shown separately) and reports associated with The transcription levels of GMFB in male normal, female normal, male HCC and female HCC tissues were analyzed by the NetworkAnalyst. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001, ns, no significance.
HCC of TOP100 DEGs in overlap, male-specific and femalespecific groups were summarized in Supplementary Table 7.

Functional Enrichment Analysis of GMFB in Patients With HCC Based on Gender
To elucidate the possible mechanism of underling GMFB expression level on the clinical outcome in HCC patients with different gender, male and female common, male-specific and female-specific GMFB co-expression genes were further analyzed. Venn diagram of GMFB co-expression genes (FDR (BH-corrected p-value) <0.01, P<0.01) in male and female HCC groups was illustrated that there were a larger number of regulated genes in males (6328 DEGs) than in females (2899 DEGs), and showing a substantial overlap (2513 DEGs) between males and females ( Figure 5A). DAVID analysis software was utilized to identified functional enrichment terms, including GO (BP, MF, and CC) and KEGG of GMFB co-expression genes in overlap ( Figures 5B, C), male-specific ( Figures 5D, E) and female-specific ( Figures 5F, G), respectively. The enrichment analysis disclosed the distinct role of GMFB co-expression genes in female and male HCC patients, despite increased expression of GMFB was observed in both male and female HCC patients.
Results also showed a high expression level of GMFB related to mitochondrial dynamics. To further capture the relationships between the terms, PPI network and MCODE component analysis of top 500 DEGs of overlap ( Figures 6A, B), malespecific ( Figures 6C, D) and female-only ( Figures 6E, F) were performed. Regardless of gender, co-expression genes of GMFB were mainly involved in oxidative phosphorylation and mitochondria function ( Figures 6A, B). In male patients, peroxisome, cell cycle and inosine monophosphate biosynthesis, PRPP + glutamine => IMP were enriched from male-specific co-expression genes of GMFB ( Figures 6C, D) while in female patients, KEGG pathways such as ribosome and ubiquitin mediated proteolysis were enriched (Figures 6E, F). Detailed information was available in Supplementary Table 8. Likewise, regardless of gender, PPI networks showed that GO annotation mainly related to MCODE 1 were mitochondrial complex and transferase complex. In males, GO annotation is mainly related to MCODE 1 of mRNA splicing, via spliceosome, RNA splicing, via transesterification reactions with bulged adenosine as nucleophile and RNA splicing, via transesterification reactions, while in females, MCODE 1 was anchoring translation, peptide biosynthetic process and ribosomal subunit. Supplementary Material reports detailed data (Supplementary Table 9).

The Increased mtDNA Copy Number and Loss of Mitochondrial Membrane Potential in Hep3B by Knocking Down GMFB
Based on the bioinformatic analysis above, mitochondria were identified to be a major target of GMFB and GMFB coexpression genes. Hence, the expression levels of mitochondrial expansion program and homeostasis-associated genes were analyzed. Results revealed that 52 genes required for mitochondrial ribosome function were down-regulated ( Table 1). Mitochondrial dynamics-associated genes (MTP18, FIS1), mitochondrial protein import-associated genes (TOMM40, TOMM7), mitophagy-associated genes (PINK1) and glucose metabolism-associated genes (FBP1, TPI1 and PCK1) were significantly down-regulated in HCC patients with high GMFB expression ( Table 1). Finally, we designed siRNA oligos against GMFB to verify the finding from the bioinformatic analysis. After the determination of the knockdown efficiency of siRNA against GMFB ( Figures 7A, B), mtDNA copy number was measured by qPCR. We found that the expression level of GMFB negatively co-expression to mtDNA copy number in Hep3B cells (Figures 7C-E). Furthermore, we also detected mitochondrial membrane potential by JC-1 assay. As shown in Figure 7F, the transition of red to green fluorescent signal was observed in GMFB knockdown HCC cells, which indicated the loss of mitochondrial membrane potential is associated with low expression levels of GMFB in Hep3B cells.

The Inhibition of Hep3B Cell Proliferation, Migration and Invasion by GMFB Knockdown
To further verify the relationship between GMFB and the malignant behavior of the HCC, we also determined cell proliferation, invasive and migratory abilities by woundhealing and Transwell assays. Post-transfection 48 hours of siGMFB oligos in Hep3B cells, we found that decreased expression of GMFB inhibited wound closure significantly ( Figures 7G, H). Transwell migration assay showed less motile with crystal violet staining in GMFB knockdown Hep3B cells ( Figures 7I, J). Transwell-matrigel invasion assays showed that GMFB knockdown inhibited Hep3B invasion ( Figures 7I, K). As shown in Table 2

DISCUSSION
In the present study, we performed a comprehensive bioinformatic analysis on the association between GMFB and HCC. We found that GMFB was highly expressed in several types of cancer including HCC ( Figures 1B, C, 2A), and confirmed the expression level of GMFB was closely associated with the TNM stage, clinical stage, and poor survival rates in HCC (Figures 1D, E). GMFB was considered to be a novel potential biomarker and therapeutic target HCC (Figure 3). Moreover, knockdown of GMFB with siRNA effectively inhibited the proliferation, migration and invasion in Hep3B cells and down-regulated the expression level of some of the MMPs. Targeting GMFB may represent a promising therapeutic strategy for HCC. However, the potential mechanisms underlying GMFB contributing to HCC remain unclear.
To clarify the potential mechanism, we found a large number of DEGs were co-expressed with the high level of GMFB, and most of them were highly associated with mitochondrial functions ( Figure 4D). Thus, we compared transcriptional profiles of DEGs related to mitochondrial protein synthesis (MRPLs and, MRPSs), Mitochondrial dynamics (MTP18 and FIS1), Mitochondrial protein import (TOMM40 and TOMM7), mitophagy (PINK1) and glucose metabolism (FBP1, TPI1 and PCK1) (44). The result indicated that a high level of GMFB down-regulated 53 genes required for mitochondrial ribosome function ( Table 1), suggesting that the high level of GMFB may suppress mtDNA replication, mitophagy, and energy metabolism in HCC. Increased expression of GMFB rather than AFP was an independent risk factor for HCC ( Figure 3E). Moreover, knockdown of GMFB with siRNA effectively inhibited the proliferation, migration and invasion in Hep3B cells and down-regulated the expression level of some of the MMPs.
The mtDNA copy number gained more and more attention in cancer research. Previous studies showed decreased copy numbers in cancer samples in HCC, bladder cancer, breast cancer, kidney clear cell carcinoma and myeloproliferative neoplasm, and increased copy number detected in chronic lymphocytic leukemia, lung squamous cell carcinoma and pancreatic adenocarcinoma (45)(46)(47)(48)(49)(50). Furthermore, the decrease in mtDNA copy number was more significant in female HCC patients than in males (51). In our study, we found that the decreased expression of GMFB leads to up-regulation of mtDNA copy number, which means the higher level of GMFB negatively regulates the amount of mtDNA copy number (Figures 3C-E). Moreover, our result also suggested that low expression levels of GMFB induced loss of mitochondrial membrane potential in Hep3b cells ( Figure 7F). Previous works showed that  mitochondria is a central regulator of the decision between cellular survival and demise (52). Mitochondrial membrane potential is an indicator of mitochondrial function, and loss of mitochondrial membrane potential resulting in the intrinsic apoptotic pathway (53). To date, there is no study has explored the role of GMFB in mitochondrial dysfunction. which might partially contribute to the malignancy of HCC. Therefore, these results suggested that GMFB may regulate mtDNA function and consequently be involved in HCC tumorigenesis and progression.  It is well known that cancer cells are commonly exposed to nutritional deficiency and hypoxia, and these factors adversely impact cancer cell metastasis (54). Cancer cells use aerobic glycolysis to support proliferation instead of mitochondrial oxidative phosphorylation (55). Therefore, inhibition glycolysis in cancer cells represents a kind of therapeutic strategy (56). In our study, GMFB down-regulated the expression level of the negative regulator of glycolysis, fructose-1,6-bisphosphatase (FBP1), Triosephosphate isomerase (TPI1) and phosphoenolpyruvate carboxykinase 1 (PCK1) ( Table 1). This suggests that high expression levels of GMFB co-expression with the promotion of glycolysis and suppression of gluconeogenesis. Furthermore, previous works showed that suppressed FBP1, a negative regulator of aerobic glycolysis, led to promoting HCC growth and metastasis (57). Knockout PCK1, a step limiting enzyme of gluconeogenesis, markedly enhanced the global O-GlcNAcylation levels, which is an emerging hallmark of HCC (58). TPI1 functioned as a tumor suppressor in HCC (59). In conjunction with down-regulated KEGG pathways from co-expression genes of GMFB, we favored that targeting GMFB/mtDNA/glycolysis may represent a novel therapeutic strategy.
According to epidemiological and clinical characteristics of HCC, there is a remarkable gender disparity of HCC incidence, primarily males. Liver cancer ranks fifth in terms of global cases and second in terms of deaths for males. Women appeared to have better survival rates than men in HCC (60). In 2018, the distribution of cases in both sexes, males, and females were 6th, 5th and 9th, and for deaths is 4th, 2nd and 6th respectively (1). When we did KM Plotter survival analysis, we found that high GMFB expression was significantly associated with poor OS of HCC patients (P=0.0042) ( Figure 3A). Unexpectedly, there was a significant OS difference between male patients with low and high expression of GMFB (P=0.00014) ( Figure 3B). However, it did not occur in females with HCC (P=0.36) ( Figure 3C).
It is widely recognized that the liver is one of the main responsible organs for estrogens (61,62). Previous studies indicate that estrogen and estrogen receptors played a role in the development and progression of HCC (62)(63)(64). But, the role of estrogen and estrogen receptors of HCC is disputed (65). It seems that increased estrogen synthesis and variants estrogen receptor in liver leads to an increased risk of HCC (66,67). On the contrary, estrogen and estrogen receptors can decrease the malignancy of HCC by arresting cell cycle progression and promoting apoptosis (68). In addition, distinct gender differences were observed in our study. GO and KEGG pathways enrichment analyses of DEGs (Figures 5, 6). Further analysis showed that more DEGs co-expressed with GMFB were regulated in the males than in the females ( Figure 5A). Results of survival analysis and current research of top 100 DEGs were shown in Supplementary Table 7. These distinct results above indicated that sex hormones might play a key role in regulating the gene expression profiling of co-expression genes of GMFB. Among top 100DEGs, some of them with statistical significance of survival analysis, such as Dodecenoyl-CoA isomerase (DCI, p value=4.3e-5 in males), G2/M phase-specific E3 ubiquitinprotein ligase (G2E3, p value=9.5e-6 in males), LSM domaincontaining protein 1 (LSMD1, p value=0.0001 in males), Neural precursor cell expressed developmentally down-regulated protein 1 (NEDD1, p value=0.0001 in males) and Ras-related GTP-binding protein C (RRAGC, p value=0.0001 in females) have not been reported in HCC yet. These genes merit further investigation.
It is well known that the MMPs play a vital role in cancer invasion and metastasis (69). We also found that expression patterns of co-expression MMPs genes of GMFB showed a gender disparity ( Table 2), which may explain the significant correlation between OS and GMFB in male HCC, not in females. Knockdown of GMFB effectively inhibited the expression levels of MMPs. Our results supported that the expression level of GMFB in male HCC modulated a subset of MMPs expression, finally contributing to male OS.
Certain limitations exist in our study. The first limitation is that we lack validation of our novel findings with human HCC samples. Since we could not obtain the levels of sex hormones, the expression level of sex hormone receptors in HCC patients used in this study, the relationship between levels of sex hormones, sex hormones reporters and GMFB expression merits further investigation. The second limitation is that mechanistic study is preliminary and should be deeply explored. In summary, this work will shed deep light onto GMFB function in HCC. Targeting GMFB may represent a promising therapeutic strategy for HCC patients. GMFB may be an invaluable diagnostic and prognostic biomarker for HCC from the perspective of precision medicine.

AUTHOR CONTRIBUTIONS
Among the authors in the list, WS designed the study, performed experiments and wrote the manuscript. CH performed part of data analysis and interpreted the experimental data. TW performed part of data analysis. JW, JPZ, FG, QO, HT, CJ, JX, and JFZ discussed the results. LL and G-TX encouraged us to investigate and supervised the findings of this work. All authors contributed to the article and approved the submitted version.

FUNDING
This work was supported by grants obtained from the Ministry of Science and Technology of China (2017YFA0104100), and the National Natural Science Foundation (81670867, 81372071, 81770942).

ACKNOWLEDGMENTS
We thank all members of our laboratory for many helpful discussions.