Hai-Hang Zhang, Jiang-Zheng Zeng, Xu Zhu, Hua-Mao Sun, Lu Yang, Yu Fu, Yan-Da Lu?

1. Department of Oncology, First Affiliated Hospital of Hainan Medical College, Haikou 570102, China

2. School of Biomedical Information and Engineering, Hainan Medical College, Haikou 571199, China

Keywords:Hepatocellular carcinoma Epigenetic factor microRNA TCGA

ABSTRACT Objective: To explore the targeting relationship between miRNA and epigenetic factors related to the prognosis of hepatocellular carcinoma (HCC), and to identify the potential impact of miRNA targeted epigenetic factors on the prognosis of HCC. Methods: The mRNA and miRNA sequencing data of all HCC samples were downloaded from the tumor Genome Atlas(TCGA) database, and the R program was used to analyze the difference of the sequencing data. The data of survival time, survival status and differential expression of miRNA were combined, and the risk score model of miRNA was constructed by univariate and multivariate Cox regression. The miRNA target genes and all the corresponding Epigenetic factors were predicted, and the differentially expressed Epigenetic factors (DEEFs) were screened.Then,the regulatory network of miRNA targeting deefs was established, the apparent factors in the network were enriched and analyzed, and the core genes in the protein-protein interaction(PPI) network and selection network were constructed. Finally, the relationship between the apparent factors and the prognosis of patients was analyzed and verified by Kaplan Meier (K-M)method. Results: 305 differentially expressed miRNAs were identified using EDGE algorithm.After Cox analysis, hsa-miR-139-5p, hsa-miR-101-3p and hsa-miR-7-5p (miR-139-5p, miR-101-3p and miR-7-5p) were finally screened as the Overall survival of HCC patients (Overall survival,OS). In addition, 34 DEEFs targeted by miRNAs were identified, among which EZH2, PKM, HJURP and CHEK1 had a significant impact on the survival of hepatocellular carcinoma. Conclusion: In this study, we successfully established a prognostic model of hepatocellular carcinoma miRNA-targeted epigenetics, and screened out epigenetic factors that are significantly related to the prognosis of hepatocellular carcinoma. It provides new potential prognostic biomarkers and therapeutic targets for HCC treatment, and lays a theoretical foundation for the follow-up basic research of HCC.

1. Introduction

Hepatocellular carcinoma (HCC) is one of the most common malignant tumors in the world and the main cause of death from cancer[1,2]. In China, chronic hepatitis and cirrhosis caused by hepatitis B virus (HBV) or hepatitis C virus (HCV) infection are the main risk factors for HCC[3,4]. Although various clinical treatments are continuously used for the treatment of liver cancer, most HCC patients have a high recurrence rate and poor prognosis, with an overall survival (OS) of only 3%-5%[5]. On the one hand, liver cancer is highly heterogeneous; on the other hand, the biological process of the occurrence and development of HCC is very complex. So far, alpha-fetoprotein (AFP) is the most widely used biomarker for liver cancer. In addition, due to the poor reliability of AFP, finding new molecular markers related to the prognosis of HCC is of great significance for the diagnosis, prognosis and treatment of liver cancer[6,7].

The occurrence and development of HCC are not only related to genetic mutations encoded by DNA, but also related to the abnormal accumulation of epigenetic regulatory genes, the activation of tumor genes, and the inactivation and loss of tumor suppressor genes[8]. Contrary to gene change, epigenetic regulation is a kind of gene expression and regulaton that does not involve changes in DNA sequence. Examples include post-transcriptional regulation of miRNA, covalent histone modification, and DNA methylation to influence chromatin structure and transcription of genes[9-10]. Although genetic changes cannot be reversed, epigenetic modification is a dynamic and reversible process[11]. The current question to be addressed by the HCC Institute is the mechanism between the heterogeneity and complexity of the genetic and epigenetic changes encoded by these DNAs and their impact on patient outcomes. Analyzing the mechanisms of complex genetic and epigenetic changes and the development of HCC may help to find effective treatments.

Although the study of epigenetics has made great progress in recent years, various drugs targeting epigenetic regulatory proteins have been developed, which can restore malignant cancer cells to a normal state[12]. However, most studies focus on a single epigenetic aspect, and there is still a lack of research on evaluating the relationship between epigenetic genes post-transcriptionally regulated by miRNAs in the prognosis of HCC. Therefore, exploring the regulatory relationship of epigenetic factors targeted by miRNAs in HCC and identifying epigenetic factors that are significantly related to prognosis can provide new ideas for the basic research of epigenetics in HCC.

2. Materials and methods

2.1 Screening of Differential Genes in HCC

Download RNA sequencing data of 371 HCC samples and 50 normal liver tissue samples and miRNA sequencing of 372 HCC samples and 50 normal liver tissue samples from the TCGA database(https://portal.gdc.cancer.gov/) . The gene expression levels of HCC samples and normal liver tissue samples were compared using R procedure. And according to the criteria of ∣log 2 FC∣> 1 and FDR <0.05, the differentially expressed miRNA and mRNA were screened.

2.2 Construction of miRNA-DEEFs regulatory network

Clinical data of all samples were downloaded from TCGA database, and the collected clinical data were integrated into a matrix with differential analysis miRNA. Univariate Cox regression analysis was used to obtain miRNAs with significant prognostic differences,and multivariate Cox regression analysis adopted stepwise regression. In addition, all epigenetic factors were downloaded from the epigenetic factor database and then interacted with miRNA target genes to identify the epigenetic factors targeted by miRNA in the model. The target genes of miRNA in the model were predicted by online database. The epigenetic factors targeted by miRNAs and differentially expressed mRNAs in HCC are intersected to screen out DEEFs. According to the regulatory mechanism of miRNA on target genes, the regulatory network of miRNA and DEEFs was finally successfully constructed.

2.3 Enrichment analysis of DEEFs and construction of protein interaction network

Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes(KEGG) pathway analysis were performed on DEEFS in the regulatory network by using the software package"Clusterprofiler" (version 3.14.3)[13]. GO analysis revealed the role of DEEFs in biological processes, cell composition and molecular functions, and KEGG analysis showed the enrichment pathway of DEEFs. STRING tool was used to construct protein-protein interaction (PPI) network to explore epigenetic gene interaction[14].In addition, core genes were selected using Cytoscape software(version 3.7.1)[15] and CytoHubba plug-in.

2.4 Comprehensive analysis of risk scoring model and clinical characteristics

In order to further explore the correlation between miRNA and the clinical characteristics of HCC in the prediction model, the Univariate COX proportional hazard regression analysis was used to identify the key clinical indicators that affect the OS rate.Clinical indicators with P < 0.05 were selected and multivariate Cox proportional hazard regression analysis was performed to determine whether the risk score model could independently predict prognosis.Draw forest plots of Univariate Cox regression analysis and multivariate Cox regression analysis to study the value of clinical indicators and the potential application of risk scoring models in clinical practice.

2.5 Survival analysis of DEEFs in the regulatory network

Firstly, the downloaded clinical data were collated, and the cases with no OS or OS less than 30 days were deleted from the HCC clinical information. Then, the patient's survival time and survival status were combined with the differentially expressed mRNA sequencing data in HCC. The KaplanMeier (K-M) method was used to analyze the relationship between the epigenetic factors in the network and the prognosis of HCC.

2.6 Statistical processing

All data were analyzed and processed by R software package(version 3.6.3) and Perl script tool. The t-test and analysis of variance (ANOVA) were used for continuous variables, and Pearson's chi-square test and Fisher's exact test were used for classification variables to detect statistical differences. When P<0.05,the difference is statistically significant.

3. Results

3.1 Difference analysis of HCC sequencing data

RNA-seq data of HCC were obtained from the Tumor Genome Atlas (TCGA) database with download time of 2020-10-18 in 374 cases of primary HCC (3 cases of recurrent tumors), of which 371 cases were untreated HCC samples and 50 normal tissues.We also downloaded the miRNA profile data of 375 cases of primary hepatocellular carcinoma (3 cases of recurrent tumors), including 372 newly treated HCC samples and 50 normal tissue samples.All tumor samples were derived from before treatment. Genes are annotated with the gencode. v22. annotation. gtf. gz probe, and if there are multiple probes corresponding to the same gene, the average value is taken. The normalized form of gene expression is log2 (TPM+1). All data involved in this research are downloaded from TCGA, and data collection and application are carried out in accordance with TCGA's published guidelines and data access strategies, and no ethical review is required. The difference analysis was performed using the "edge" R software package [16] in BioConductor on R (version 3.6.3). According to (LogFC>=1, FDR<0.05), the sequencing data of differentially expressed mRNA was screened, and 25 differentially expressed mRNAs with high and low expression were selected to draw a heatmap by cluster analysis (Fig.1A). The difference between groups of 4815 up-regulated genes and 1367 down-regulated genes in mRNA was shown by drawing a Volcanic plot (Fig. 2A). The same method was used for differential analysis of miRNA expression profile data, and 25 differentially expressed miRNAs with high and low expression were selected to draw a heatmap through cluster analysis (Fig. 1B), of which 257 genes were up-regulated, 48 genes were down-regulated, and a Volcanic plot was drawn (Fig. 2B) .

Figure1 Heatmap of differentially expressed mRNA and miRNA in HCC(Column: representative sample; row: mRNA or miRNA; red: high expression; green: low expression)

Figure2 Volcanic plots of differentially mRNA and miRNA in HCC(Gray:no differentially expressed genes; red: up-regulated differential genes; blue:down-regulated differential genes)

3.2 Construction and selection of predictive models

After obtaining candidate miRNAs by the above method, the entire TCGA-HCC data set (n=373) was divided into a training set (172 cases) and a test set (171 cases) using random numbers generated by the computer. Then the clinical data of HCC patients were sorted out, and the cases with no OS or OS less than 30 days were deleted, and then the patient survival time, survival status and differential miRNA expression data were merged. The model is built in the training set, and then applied to the test set and the entire TCGA-HCC data set to select the best model. Using univariate Cox analysis, miRNAs with prognostic value for overall survival(OS) of HCC were selected according to P<0.01. Multivariate Cox regression analysis established a good and stable prediction model based on 3 miRNAs in the training set（Table.1). The Kaplan-Meier (K-M) survival curve indicates that the HR value of miR-7-5p is greater than 1, which is a high-risk miRNA.With the increase of miR-7-5p expression, the prognostic risk of patients gradually increased (Fig. 3C). On the contrary, miR-139-5p and miR-101-3p have HR values less than 1, which are low-risk miRNAs. With the increase of their expression levels, the prognosis risk of patients decreased gradually(Fig. 3A, B). The stepwise weighted correlation coefficients of multivariate Cox regression are as follows: Survival risk score = (-0.29371 × ExpmiR-139-5p) + (-0.24725 × ExpmiR-101-3) + (0.376329 × ExpmiR-7-5p). According to the median of the risk scores of patients with HCC, the patients in the training group were divided into two groups with high and low expression.The OS of patients in the high-risk group (n=86) was worse than that in the low-risk group (n=86) (Fig. 4A). We use time-related ROC to evaluate the reliability of the risk scoring model, and the area under the curve (AUC value) is 0.742 (Fig. 5A), indicating that the prediction model has achieved good accuracy in survival monitoring.At the same time, the risk score curve, survival status and expression pattern showed that patients with lower scores in the two groups usually had a better prognosis than patients with higher risk scores(Fig 6A). The accuracy of the prediction model is further proved.

Figure3 Kaplan-Meier survival curves of differentially expressed miRNA associated with prognosis in HCC(A: Relationship between miR- 139-5p expression in HCC tissue and prognosis; B: Relationship between miR-101-3p expression in HCC tissue and prognosis; B: Relationship between miR-7-5p expression in HCC tissue and prognosis)

Table1 List of miRNAs in HCC prognostic models

3.3 Evaluate the predictive model

To further confirm the risk scoring model, it was verified on the test set and the entire TCGA-HCC data set. The model successfully divided HCC patients in the test set (n=171, P=6.561e-03) into a high-risk group (n=90) and a low-risk group (n=81), and there were significant differences in OS (Figure 4B). ROC analysis over time showed that the model performed well in the prediction of HCC patients in the test group. The RUC value for this group was 0.757 (Figure 5B). Similarly, validation on the entire TCA-HCC dataset (n=343, P=8.273e-06) showed that patients in the high-risk group (n=176) had poorer OS compared to patients in the low-risk group (n=167).(Figure 4C). ROC analysis over time shows that the model also has a good predictive performance for HCC patients in the entire TCGA-HCC data set. The RUC value for this group was 0.757 (Figure 5C). The distribution of risk scores, survival status,and expression patterns in the test set and the entire TCGA-HCC data set also showed consistent results with the training group. The survival rate of patients with higher risk scores was lower than that of patients with lower risk scores (Fig. 6B, C). The results not only prove that there are significant differences in survival time and risk degree between high and low risk HCC patients, but also prove the accuracy of the risk model.

Figure4 Kaplan-Meier survival curve of miRNA risk model related to prognosis in HCC

Figure5 Receiver operating characteristic (ROC) curve analysis on the Accuracy of miRNA Risk Models in HCC to predict the prognosis of patients

Figure6 Risk score, survival status and gene expression profile heat map of HCC patients

3.4 Comprehensive analysis of risk scoring model and clinical characteristics

To assess whether the predicted miRNA (miR-139-5p, miR-101-3p, and miR-7-5p) interactions are related to the overall survival of HCC patients. First, according to the miRNA risk score, samples of HCC patients are divided into high-expression groups and lowexpression groups. Then, determine whether the prognostic value of the predictive signal is independent of other clinicopathological variables in HCC patients. Our selected variables included age,sex, TNM stage of pathology, pathology grade, and AJCC stage,and our risk score for univariate and multivariate Cox regression analyses. The results indicate that the predictive risk score model is an independent prognostic factor for HCC patients (Fig. 7). In summary, the predictive risk score predicts the prognosis of HCC patients independent of other clinical characteristics.

Figure7 Evaluation of the prognostic accuracy of miRNA risk models and patients' clinicopathological variables in HCC (A; Cox regression analysis of univariate clinical phenotype; B: Cox regression analysis of multivariate clinical phenotype)

3.5 Construction of miRNA-DEEFs regulatory network in HCC

Using miRTarBase, TargetScan and miRDB databases to predict miRNA target genes in the model, a total of 21763 pairs of miRNAtarget genes were obtained. First, 720 epigenetic factors (EpiFactors database (http://epifactors.autosome.ru/) and miRNA target genes were intersected from the epifactor database, and 533 epigenetic factors targeted by miRNA in the model were identified. Then, the epigenetic factors targeted by miRNAs in HCC and the differentially expressed mRNAs were crossed to screen out 63 DEEFs. Among them, 34 epigenetic factors were negatively correlated with miRNA expression (8 down-regulated target genes, 26 up-regulated target genes). Finally, based on the 3 miRNAs and 34 differentially expressed epigenetic factors related to the prognosis of HCC, a regulatory network of miRNA-DEEFs was constructed (Fig. 8).

Figure8 Screening and analysis to obtain the miRNA-DEEFs network diagram in HCC (Triangle; representative miRNA; ellipse: epigenetic factor;red: high expression; green: low expression)

3.6 Functional enrichment analysis of DEEFs and protein interaction network

GO enrichment analysis is performed on DEEFs targeted by miRNAs. GO enrichment consists of three parts: Biological process(BP), Cellular component (CC) and Molecular function (MF). The first three GO terms of biological processes are covalent chromatin modification, histone modification and peptidyl lysine modification.The first three GO terms of cell components are chromosome region, nuclear chromatin and chromosome, centromere region, The first three GO terms for molecular functions are methyltransferase activity, transferase activity, transfer of one-carbon group and S-adenosylmethionine-dependent methyltransferase activity (Fig.9A). KEGG enrichment results confirmed that lysine degradation and MicroRNAs in cancer are the most significant enrichment pathways (Figure 9B). The enrichment of GO and KEGG showed that the predicted target genes are closely related to the regulation of epigenetic factors. It mainly affects covalent modification of chromatin, histone modification and peptidyl lysine modification,methyltransferase activity, lysine degradation and regulation after miRNA expression.In order to further explore the interactions between epigenetic factors, we use the search tool (STRING, https://string-db.org/) of protein-protein interaction network (PPI) analysis to construct a miRNA-targeted table View factor PPI network(Fig.10A), use Cytoscape software (version 3.7.1) and CytoHubba plugin to select core genes for visualization(Fig. 10B).

Figure9 Enrichment analysis of miRNA regulated epigenetic factors in HCC(A:GO;B:KEGG)

Figure10 Regulatory network of differential apparent factors (A: PPI network; B: network core module)

3.7 Survival analysis of DEEFs regulated by miRNA

KaplanMeier (K-M) method was used to analyze the relationship between the expression level of epigenetic factors regulated by miRNA and the prognosis of HCC. When the expression levels of miR-139-5p and miR-101-3p decrease in HCC, the high expression of their regulated epigenetic factors (EZH2, PKM, HJURP and CHEK1) is associated with a poor prognosis. It suggests that these epigenetic factors may play a role in cancer (Fig.11A-D). Although it has been confirmed above that the high expression of miR-7-5p indicates a poor prognosis in HCC patients, in the K-M survival analysis of this study, no antitumor effect of the epigenetic factors targeted by miR-7-5p on the prognosis of HCC patients was found.

Figure11 Kaplan Meier survival curve of differential apparent factors in HCC (A: Compared with the low expression group, the overall survival rate of liver cancer patients with high expression of EZH2, PKM, HJURP and CHEK1 was significantly reduced)

4. Discussion

The expression characteristics of miRNA in tumors are different from those in normal tissues, and also vary with tumor types.Increased miRNA expression is often observed in HCC, and these up-regulated miRNAs act as oncogenic miRNAs in HCC[17].Therefore, when the regulation of miRNA is released, it can affect the proliferation, invasion and metastasis of tumor cells, which may be the target of tumor suppressor genes[18]. Obviously, the change of miRNA endogenous expression is an important mechanism of liver cancer. Considering the synergistic effect of miRNA on epigenetic regulation after transcription, we tried to study the influence of miRNA and targeted epigenetics in HCC through integration analysis. Finally, a miRNA-targeted epigenetic regulation network was constructed to comprehensively analyze the synergistic effects of miRNA on epigenetic regulation after transcription, and the results showed that miR-139-5p, miR-101-3p and miR-7-5p were significantly correlated with the prognosis of HCC. In addition,survival analysis of miRNA-targeted epigenetics found that EZH2,PKM, HJURP, and CHEK1 have a significant impact on the survival of HCC.

Studies have found that miR-139-5p has anti-cancer effects in different organs, and high expression of miR-139-5p is associated with a better prognosis[19]. Through the analysis of HCC sequencing data, we found that miR-139-5p in liver cancer tissues was a significantly up-regulated miRNA compared to normal tissues.However, the overexpression of miR-139-5p can reduce the invasion and proliferation ability of HCC cells by increasing the expression of the SLATRK4 pathway [20], suggesting that miR-139-5p will become an effective target for the treatment of HCC.There is evidence that the high expression of the epigenetic factor HJURP regulated by miR-139-5p promotes the proliferation of HCC cells. Clinically, the high expression of HJURP is related to the unfavorable prognosis of individuals with HCC[21]. The role of PKM in liver cancer is to enhance glycolysis and promote the progress of HCC. The main mechanism is to inhibit the reprogramming of HCC glucose metabolism, cell proliferation and metastasis[22-23]. Therefore, high expression of PKM promotes progression in HCC. In this study, it was confirmed that miR-139-5p and targeted regulatory epigenetic factors HJURP and PKM play an important role in the occurrence and development of human HCC. The newly identified epigenetic factors HJURP and PKM may become new targets for the treatment of HCC.

In addition, hepatitis B virus (HBV) can inhibit the expression of miR-101-3p by inhibiting the activity of its promoter. In HCC,the expression of miR-101-3p is significantly lower than that of normal tissues[24]. The specific miR-101-3p inhibitor can enhance the cell proliferation, metastasis and invasion ability of HCC[25-26]. The epifactor EZH2 targeted by miR-101-3p is one of the most significant deregulation in HCC, and the up-regulation of EZH2 is closely related to the progression and metastasis of HCC[27]. In addition, CHEK1 overexpression in hepatocellular carcinoma leads to poor overall survival (OS) in HCC patients[28]. In this study, the expression of miRNA-101-3p in HCC was up-regulated and the patient's prognosis was poor. In addition, patients with upregulated epigenetic factors EZH2 and CHEK1 have a poor prognosis,confirming that miR-101-3p and its targeted regulation of epigenetic factors EZH2 and CHEK1 play an important role in the etiology of human HCC.

The discovery of miR-7-5p in this study expands our understanding of post-transcriptional epigenetic factor regulation during cancer development. It is well known that miR-7-5p can be directly bound by the long non-coding RNA (lncRNA) RUSC1-AS1 to promote the proliferation and reduce apoptosis of HCC cells.The upregulated expression in hepatocytes can promote the occurrence and development of HCC in vivo [29].This suggests that miR-7-5p plays a critical role in the development of HCC. No significant correlation was found between miR-7-5p targeting epigenetic factors and HCC prognosis. However, we speculate that they play an important role in the occurrence and development of HCC.

Although the importance of miRNAs is recognized in regulating the expression of protein-coding genes, the precise functions of miRNAs are still difficult to determine. However, miRNAs play an important role in human carcinogens. A large number of studies have found that miRNAs can act synergistically on target genes.Focusing on the interaction between miRNAs and target genes may more truly capture the underlying pathological mechanism of liver cancer. In this study, the successfully constructed miRNA targets the regulatory network of DEEFs, revealing the relationship between the differentially expressed miRNAs in HCC cancer and the prognosis of HCC. Among them, miR-139-5p and miR-101-3p, which are related to tumor inhibitors, are reduced in expression during tumorigenesis, development and metastasis, and may be potential therapeutic targets. On the contrary, the high expression of miR-7-5p, which is related to cancer genes, may promote the progression of HCC. The regulation of various biological processes by miRNAs in normal cells is strictly monitored[30]. Each miRNA may independently control epigenetic regulation through various regulators[31-33]. Therefore, the construction of a regulatory network of miRNAs targeting differential epigenetics reveals the complex regulatory mechanism of hepatocellular carcinoma epigenetics, and helps to deepen the study of miRNAs in HCC on the regulation of epigenetic genes.

Conflict of interest

All authors declare that there is no conflict of interest