Introduction
Liver cancer is currently the sixth most common malignant tumor with the fourth highest mortality rate worldwide. Hepatocellular carcinoma (HCC) accounts for about 80% of liver cancer and its incidence is still increasing worldwide.1 Even though the current treatment system based on surgical resection has been developed, the 5-year overall survival remains relatively poor, because of the difficulties associated with early diagnosis and metastasis.2–4 Therefore, it is essential to elucidate the underlying molecular mechanisms for the development and progression of HCC and, in particular, to identify potential molecular biomarkers for the prediction of the prognosis of patients with HCC. At present, alpha-fetoprotein, glypican-3, and other proteins have been evaluated, but the sensitivity and specificity are low because of the heterogeneity of HCC.5
Along with the rapid development of mRNA-sequencing technologies and bioinformatics data, a novel approach integrating gene signatures with clinical parameters has been developed with a great advantage in cancer prognosis.6–9 Among the gene signatures, those based on the immune-related genes are emerging as a novel hallmark.10 It has been demonstrated that HCC is an immunogenic tumor, where the tumor immune microenvironment is one of the factors that influences tumor initiation and progression.11–13 Accumulating evidence shows that the tumor immune system in HCC is involved in the progression and prognosis. Thus, a prognostic model that is derived from immune-related-genes is expected to be a powerful tool to predict prognosis, and will help physicians to better monitor patient treatment outcomes and improve the overall survival rate. Indeed, immune-related genes have been used to establish prognostic models for various types of cancers, including cervical, ovarian, and lung cancer, with good performance,14–16 and which support the potential of immune-related genes to become tumor prognostic markers.
Nuclear receptor subfamily 6, group A, member 1 (NR6A1), an orphan member of the nuclear receptor superfamily, is known to regulate the cell cycle and to affect epithelial-mesenchymal transformation in cancer.17 At present, available data on the prognostic and functional roles of NR6A1 are mainly related to reproductive diseases. Sequence comparisons indicate that NR6A1 is highly conserved among species.18 Studies of mouse embryonic development have shown that NR6A1 is involved in several distinct developmental processes, including possible primordial germ cell differentiation and early development, late development of the gastrula, body axis formation, and neurogenesis.19NR6A1 knockdown increases lipid accumulation and insulin-induced proliferation and migration of HepG2 cells.20 However, few studies on the prognostic roles of NR6A1 in cancer have been performed. It has been reported that hypomethylation of NR6A1 is associated with poor overall survival of acute myeloid leukemia in the Southwest Oncology Group (SWOG) and Cancer Genome Atlas (TCGA) cohorts,21 and that overexpression of NR6A1 restores the proliferation, migration and invasion in gastric cancer cells lacking hsa_circ_001653.22 Recently, a five-gene (HDAC1, BIRC5, SPP1, STC2, and NR6A1) prognostic model was established for the prediction of overall survival in patients with HCC.23 However, the performance of the model was only of acceptable predictive ability, and there was no in vitro or in vivo validation.23 Moreover, the expression profile and molecular biological functions of NR6A1 during the development and progression of HCC, and its prognostic significance remain unknown.
Therefore, the aim of the present study was to establish and validate a prognostic model for HCC by identifying immune-related differentially expressed genes (IR-DEGs) through the TCGA Liver Hepatocellular Carcinoma (TCGA-LIHC) dataset. In addition, as NR6A1 is one of 10 identified genes included in the prognostic model, but is relatively less extensively investigated compared with other nine IR-DEGs,24–26 this study further investigated the potential role of NR6A1 in the progression of HCC.
Methods
Data sources
mRNA expression and clinicopathological data were downloaded from the Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) dataset (https://tcga-data.nci.nih.gov/tcga/ ). Immune-related genes were obtained from the ImmPort database (www.immport.org ), which includes immune categories such as interleukins, tumor necrosis factor family receptors, T/B-cell receptor signaling pathway, and so on. Perl (ActivePerl 5.28 https://docs.activestate.com/ ) was used for data extraction and matching. Transcriptional data of 371 HCC cases and 50 healthy volunteers from a TCGA-LIHC dataset were extracted and the matched genes were identified. The main clinical characteristics of HCC cases are summarized in Supplementary Table 1, which is available online (https://portal.gdc.cancer.gov/ ). The protocol of the present study was approved by the Institutional Review Board of Qingdao Municipal Hospital, and patient consent was waived as the data used in the study were open to the public. The study was conducted following the ethical guidelines of the Declaration of Helsinki, 2013 revision.
Identification of IR-DEGs
The Limma package of R software was used to screen IR-DEGs in HCC and normal tissues and the mRNAs with p-values <0.05 and |Log fold change (FC)| >2 were defined as IR-DEGs, and included in the subsequent analysis. The IR-DEGs were visualized with the pHeatmap R package and were further tested by univariate Cox regression analysis. Those with a p<0.001 were included in the Lasso Cox regression analysis to identify the gene signature for the establishment of the prognostic model.27
Establishment and validation of the prognostic model
The prognostic risk formula was derived from the combination of the expression levels of the IR-DEGs identified in the multivariate Cox regression model multiplied by the regression coefficients. The risk score was calculated as (Coefficient mRNA1 × expression level of mRNA1) + (Coefficient mRNA2 × expression level of mRNA2) + ⋯ + (Coefficient mRNA × expression level of mRNA).28
The median risk score was used as a cutoff to divide the 371 HCC patients into high- and low-risk groups. Kaplan-Meier analysis was performed to estimate and compare the difference in overall survival (OS) between the two groups. Receiver operating characteristic (ROC) curves were plotted and the area under the curve (AUC) predicting the 1- and 3-year survival was determined with the R statistics survival and survminer packages, respectively. To further verify the predictive capacity of the prognostic model, 231 HCC patients from the ICGC database were used for external validation.
Functional and pathway enrichment analysis of HCC patients with high- and low-risk scores
Gene ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes KEGG) pathway analyses were performed on the TCGA-LIHC dataset and visualized using the ggplot2 R package. Differences in biological functions and pathways between the high- and low-risk groups were compared. Gene Set Enrichment Analysis (GSEA, https://www.gsea-msigdb.org/gsea/index.jsp ) was also used to explore the potential pathways of NR6A1 in HCC.
Cell culture
Cells of a human HCC cell line, Huh7 were purchased from the Chinese Academy of Sciences (Beijing, China) and cultured in Dulbecco’s modified Eagle’s medium (DMEM, Biological Industries, Kibbutz Beit Haemek, Israel) supplemented with 10% fetal bovine serum (Biological Industries) in a 5% CO2 incubator at 37°C, and then used in the following experiments.29
Transfection
An NR6A1-expressing plasmid (Youbio Biological Technology, Hunan, China) was transfected into Huh7 cells by using Lipofectamine 3000 reagent (Invitrogen, California, USA) following the manufacturer’s instructions. Cells transfected with an empty vector (Youbio Biological Technology, Hunan, China) were used as controls. The final volume was 100 µL DMEM, with 3 µL Lipofectamine 3000 Reagent, 2 µL P3000 Reagent, and 1 µg plasmid in each well of a 12-well plate.
Western blotting
Cells were lysed by grinding the cell precipitate with 2 × SDS lysate, and the proteins were separated by sodium dodecyl sulfate-polyacrylamide gel electrophoresis and transferred to polyvinylidene difluoride membranes that were blocked with 5% milk. The primary antibodies were rabbit anti-human anti-GAPDH (1:3,000, Bioss, Beijing, China) and anti-NR6A1 (1:1,000, Proteintech, Illinois, USA) antibodies. The secondary antibodies were goat anti-rabbit IgG secondary antibodies (1:5,000, Jackson ImmunoResearch, Pennsylvania, USA). Chemiluminescence gel was used for image development.
Cell proliferation assay
Cell proliferation was determined with a cell counting kit-8 (CCK8) assay (Vazyme, Nanjing, China). In brief, cells seeded in 96-well plates were transfected with the empty vector or the NR6A1-expressing plasmid. CCK8 was then added to each well at 0, 24, 48, and 72 h. After incubation in the dark at 37°C for 1 h, absorbance was measured at 450 nm.
Colony formation assay
Cells transfected with the empty vector or NR6A1-expressing plasmid were cultured at 37°C for 24 h. After digestion with trypsin, single cell suspensions were prepared. After cell counting, 1,000 cells were seeded into each well of a 6-well plate. After 1 week in culture, paraformaldehyde and 0.1% crystal violet were added to fix and stain the colonies that had formed.
Wound healing assay
Cells were seeded into 12-well plates and transfected with the empty vector or NR6A1-expressing plasmid. After 24 h culture at 37°C, a pipette tip was used to wound in a single layer of cells. The cells were cultured in serum-free DMEM at 37°C and photographed at 0 and 24 h. The width of wound was measured and cell mobility was calculated.
Transwell migration assays
Cells transfected with the empty vector or NR6A1-expressing plasmid were digested with trypsin after 24 h culture at 37°C. The digested cells were then suspended in serum-free DMEM and cultured in the upper compartment of 24-well Transwell chambers with an 8.0 µm pore size polycarbonate membrane (Corning Incorporated, New York, USA). DMEM containing 25% fetal bovine serum was added in the lower compartment of the chambers,. After culture at 37°C for 24 h, the cells in the upper chamber were wiped with cotton swabs, and the cells penetrating the membrane were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet. The invasive cells in five random fields were counted by light microscopy.
Subcutaneous xenograft experiment
Six 6-week-old male specific-pathogen-free BALB/c nude mice (Jinan Pengyue, Jinan, China) were randomly divided into two groups. One group was injected subcutaneously in the right flank with 2×106 Huh7 cells transfected with the empty vector resuspended in DMEM, and another with 2×106 Huh7 cells transfected with NR6A1-expressing plasmid resuspended in DMEM. The mice were sacrificed after approximately 4 weeks, tumor formation was observed, and tumors, if any, were removed and weighed. The volume was calculated as length × width2/2. Animal procedures were approved by the Animal Care and Use Committee of Qingdao Municipal Hospital.
Statistical analysis
GraphPad Prism (www.graphpad.com ) and R software v3.6.3 (www.r-project.org ) were used for statistical analyses. P-values <0.05 was considered statistically significant unless specified above.
Results
Selection of IR-DEGs for generating a prognostic gene signature and model
Overall, 5,482 DEGs with |Log FC| >2 and p<0.05 were observed between the 371 HCC cases and 50 normal cases, which are shown in the volcano plot (Fig. 1A) and the generated heatmap (Fig. 1B). Of these genes, 46 IR-DEGs with a p-value of <0.001 were detected by the univariate Cox regression analysis (Fig. 1C). The detailed mRNA expression levels of the genes are shown in the heatmap (Fig. 1D). In the Lasso Cox regression analysis, 16 of the 46 genes were found to be able to construct the prognostic gene signature (Fig. 2A). Finally, 10 genes were identified in the multivariate Cox regression analysis in the TCGA-LIHC cohort (Fig. 2B). The 10 genes used to construct the signature model were heat shock protein 4 (HSPA4), fatty acid binding protein 6 (FABP6), microtubule-associated protein tau (MAPT), N-myc downstream regulated gene 1 (NDRG1), adipokine prepropeptide (APLN), interleukin-17D (IL17D), luteinizing hormone β subunit (LHB), secreted phosphoprotein 1 (SPP1), glucagon-like peptide 1 (GLP1R), and NR6A1 (Fig. 2C).
Establishment and validation of the prognostic model
The prognostic model established using the 10 IR-DEGs mentioned above was formulated as: risk score = SPA4 × 0.0000483 + FABP6 × 0.0008682 + MAPT × 0.0004037 + NDRG1 × 0.0000105 + APLN × 0.0002845 + IL17D × 0.0003986 + LHB × 0.0119960 + SPP1 × 0.0000014 + GLP1R × 0.0011598 + NR6A1 × 0.0005548. The risk scores of the 371 HCC patients ranged from 0.389 to 39.818 (Fig. 3A), with a median score of 0.823. With 0.823 as a cutoff, the patients were divided into high-and low-risk groups. As shown in Fig. 3B, the OS rate was significantly lower in the high-risk group than in the low-risk group (p<0.001), and the AUCs for predicting 1- and 3-year survival were both 0.834 in this prognostic model (Fig. 3C). In addition, patients survived longer in the low-risk group than in the high-risk groups (p<0.001; Fig. 3D). The expression profiles of the 10 genes in the two groups are illustrated as a heatmap in Fig. 3E. These results indicate that the prognostic model constructed based on the 10 IR-DEGs (HSPA4, FABP6, MAPT, NDRG1, APLN, IL17D, LHB, SPP1, GLP1R, and NR6A1) effectively predicted the prognosis of the HCC patients. In the verification of the performance stability of the prognostic model in the ICGC dataset, the AUC values for predicting 1- and 3-year survival were 0.772 and 0.738, respectively.
Differences in biological function and pathways between the high- and low-risk HCC groups
There were differences in biological function such as regulation of artery morphogenesis, and acyl-COA modulation between the high- and low-risk HCC patients as shown by GO enrichment analysis (Fig. 4A). In addition, the KEGG analysis identified significant differences in pathways such as nitrogen metabolism and fatty acid degradation, indicating the potential differences in functional genes expression in the two groups (Fig. 4B).
Upregulation of NR6A1 expression in tumor tissues of HCC patients
To further investigate the potential role of NR6A1 in the progression of HCC, we assayed NR6A1 expression in the ten IR-DEGs. The expression of NR6A1 mRNA was significantly higher in the tumor tissue of HCC patients than in the normal liver tissue of healthy volunteers in the TCGA-LIHC dataset (Fig. 5A). In addition, the expression of NR6A1 in tumor tissue increased significantly with the increase of grade/stage (Fig. 5B, C). In the UALCAN cancer database, tumor tissues were divided into groups with high (n=91) and low-medium NR6A1 expression (n=274) groups. The survival rate of HCC patients with low-medium NR6A1 expression was significantly higher than that of those with high NR6A1 expression (Fig. 5D).
Enrichment of pathways related to tumor proliferation and migration in HCC patients with high NR6A1 expression
GSEA revealed that the cell cycle, mTOR, WNT, and ERBB signaling pathways that are related to cell proliferation, migration, or epithelial-to-mesenchymal transition were significantly enriched in HCC patients with high NR6A1 expression (Fig. 6A–D).
Promotion of proliferation and migration by NR6A1 in HCC cells
The mRNA expression of NR6A1 in Huh7 cells was up-regulated by transfecting the NR6A1-expressing plasmid (Fig. 7A). CCK8 assays showed that upregulating NR6A1 promoted the proliferation of Huh7 cells (Fig. 7B). Compared with the control cells transfected with the empty vector, the number of cell colonies formed by transfected Huh7 cells significantly increased (Fig. 7C). Moreover, wound healing and Transwell migration experiments showed that up-regulation of NR6A1 expression promoted Huh7 cell migration (Fig. 7D, E).
Promotion of the growth of HCC cells by NR6A1 in vivo
The volume and weight of subcutaneous tumors in mice injected with Huh7 cells transfected with the NR6A1-expressing plasmid were significantly greater than those in mice injected with the control Huh7 cells transfected with the empty vector (1.390±0.418 vs. 2.438±0.725 mm3, p=0.048 for volume, and 1.364±0.715 vs. 2.695±0.674 mg, p=0.039 for weight). Moreover, liver metastases were found in the NR6A1-overexpressing group (Fig. 8).
Discussion
This study established and validated a risk-score based prognostic model for HCC based on ten immune-related genes, HSPA4, FABP6, MAPT, NDRG1, APLN, IL17D, LHB, SPP1, GLP1R, and NR6A1, and by analysis of a TCGA-LIHC dataset of 371 cases of HCC and 50 healthy volunteers. The prognostic model produced AUC values of 0.834 for 1- and 3-year survival in the TCGA-LIHC dataset, and 0.772 and 0.738, respectively, in the ICGC verification dataset, indicating that the prognostic model had good performance in predicting the survival of HCC patients. Moreover, the expression of NR6A1 mRNA was significantly up-regulated in HCC patients, compared with that in healthy volunteers, and NR6A1 up-regulation was significantly associated with a low survival rate of HCC patients. GSEA analysis showed that signaling pathways related to proliferation and metastasis, including cell cycle, mTOR, WNT, and ERBB signaling pathways, were significantly enriched in HCC patients with NR6A1 overexpression. In addition, in vitro and in vivo experiments showed that NR6A1 promoted cell proliferation, invasiveness, migration, and tumor formation and growth.
The study has several strengths. First, the identification of the ten immune-related gene signatures is of scientific significance and has clinical implications. It not only helps elucidate the pathogenesis of HCC, but also provides insights into future research on their roles in HCC development. Most of the ten genes in our prognostic model are known to be involved in various cancers. Among the genes, NR6A1 has been reported to be involved in the regulation of embryonic stem cell differentiation, reproduction, and neuronal differentiation.19,30NR6A1 expression is also increased in prostate cancer, where overexpression reduces G0/G1 phase cell cycle arrest and promotes metastasis and invasiveness of prostate cancer cells in vitro and the growth of prostate cancer cells in vivo.31 In this study, we observed that NR6A1 overexpression was associated with the enrichment of cell cycle, mTOR, WNT, and ERBB signaling pathways in HCC patients. It is thus conceivable that NR6A1 promotes proliferation and metastasis of HCC through the cell cycle, mTOR, WNT, and ERBB signaling pathways. HSPA4 has been reported to promote the proliferation, invasion, and metastasis of cancers, and HSPA4 overexpression is associated with a poor prognosis of HCC.32,33NDRG1 has been reported to be involved in the tumorigenesis, metastasis, recurrence, and prognosis of HCC.34IL17D has been shown to be of diagnostic and prognostic values for HCC.35SPP1 overexpression is reported to be associated with the invasion, metastasis of colorectal cancer.36 It has been shown that FABP6 expression is significantly increased in patients with colorectal cancer,37 and that FABP6 is involved in tumor cell invasion, angiogenesis, and progression of glioma.38MAPT may have a potential cell apoptosis function in cancers of the central nervous system,39 and MAPT (Tau) is associated with various pathologies including various cancers.40APLN has an important role in the development of HCC,41LHB has a key role in breast tumorigenesis,42 and GLP1R is associated with an increased risk of pancreatitis and pancreatic cancer.43 Taken together, the above mentioned ten genes are not only involved in the development of HCC, but also are prognostic factors for HCC. More important, they may be promising therapeutic targets for HCC treatment.
Second, a novel and effective prognostic model for HCC established by integrating ten IR-DEGs, may be of potential clinical implications. It may have an important role in clinical practice. For example, the prognostic model can guide treatment decisions. The risk scores of the 371 HCC patients in the TCGA database ranged from 0.389 to 39.818, with a median score of 0.823, which was defined a cutoff to differentiate high-and low-risk patients. If the risk score of an HCC patient was greater than 0.823, then it is predicted a poor prognosis, which can guide physician discussions with patients and their families on the potential survival period and provide the most appropriate treatment options.
Third, this is the first report, to our best knowledge, that explores the role of NR6A1 in the progression of HCC in both in vitro and in vivo experiments. Among the 10 identified IR-DEGs that were included in the prognostic model, NR6A1 is less extensively investigated. The roles of other nine genes in HCC progression, especially HSPA4, NDRG-1, and SPP1, have been investigated in the previous studies.24–26 In this study, both in vitro and in vivo experiments demonstrated that NR6A1 promoted HCC cell proliferation, invasiveness, migration, and formation and growth of malignant tumors.
The study has a few limitations. First, the data were derived from a public source, a TCGA dataset, and their quality and reliability largely depend on the dataset. Second, the AUC for our prognostic model was about 0.834, suggesting that there is still room to further improve the performance. Third, the study mainly focused on the prognostic role of IR-DEGs and the role of NR6A1 in the progression of HCC. Thus, the complex biological process of HCC pathogenesis and the role of nine other IR-DEGs were not explored. IR-DEGs other than NR6A1 may be also responsible for the enriched pathways identified in GSEA analysis, and NR6A1 may not be the only responsible one. Therefore, subsequent large studies are required to further improve the performance of the prognostic model and especially to uncover the molecular mechanisms by which NR6A1 promotes HCC.
Conclusions
A prognostic model based on a novel ten immune-related gene signature was established, with good performance in predicting the survival of HCC patients. NR6A1 was significantly up-regulated in HCC, and NR6A1 overexpression was associated with a poor prognosis of HCC. NR6A1 promoted cell proliferation, migration, and growth of HCC, probably through the cell cycle, mTOR, WNT, and ERBB signaling pathways.
Supporting information
Supplementary Table 1
The main clinical characteristics of HCC cases.
(XLSX)
Abbreviations
- APLN:
adipokine prepropeptide
- DEG:
differentially expressed gene
- FABP6:
fatty acid binding protein 6
- GLP1:
glucagon-like peptide 1
- GO:
gene ontology
- GSEA:
gene set enrichment analysis
- HSPA4:
heat shock protein 4
- IL17D:
interleukin-17D
- KEGG:
Kyoto Encyclopedia of Genes and Genomes
- LHB:
luteinizing hormone β subunit
- LIHC:
liver hepatocellular carcinoma
- MAPT:
microtubule-associated protein tau
- NDRG1:
N-myc downstream regulated gene 1
- NR6A1:
nuclear receptor subfamily 6 group A member 1
- ROC:
receiver operating characteristic curve
- SPP1:
secreted phosphoprotein 1
- TCGA:
The Cancer Genome Atlas
- TNM:
tumor node metastasis
Declarations
Acknowledgement
The authors thank the TCGA project and other groups for providing invaluable datasets for statistical analyses. We thank Medjaden Inc. for scientific editing of this manuscript.
Data sharing statment
The datasets used in the study are available from the corresponding author on reasonable request.
Funding
None to declare.
Conflict of interest
YNX has been an editorial board member of Journal of Clinical and Translational Hepatology since 2013. The other authors have no conflict of interest related to this publication.
Authors’ contributions
Study concept and design (SYX, YNX, LKZ), acquisition of data, experiments, analysis and interpretation of data (ZHL, JZ), drafting of the manuscript (ZHL, JZ), revision of the writing and project supervision (SYX, YNX, LKZ). All authors have made a significant contribution to this study and have approved the final manuscript.