Development an Immune-Related MicroRNA Risk Index in Hepatocellular Carcinoma

Purpose Hepatocellular carcinoma (HC) has emerged as one of the most prevalent malignancies on a global scale. Recently, immunotherapy has achieved favorable effectiveness in the management of multiple cancers. However, there are limited therapeutic options for advanced HC. As the liver is a special immune organ, we intend to uncover potential and effective immunotherapeutic modalities for HC. Our study was designed to develop specific immune-related miRNAs (IRMs) for outcome assessment and individualized strategies for the management of HC. Methods The miRNA-seq and survival data of TCGA-LIHC dataset was enrolled into this program. We first collected IRMs from Immune-miR website. Differentially expression analysis was applied to screen aberrantly expressed IRMs. In order to set up an IRM-related index (IRMRI) in HC, we conducted the Cox relevant methods. Next, the statistical approaches (survival curve and ROC curve analyses) were utilized to detect the evaluation capacity of our IRMRI. Subsequently, we obtained the target genes of hub miRNAs from IRMRI through three miRNA-related predictive online tools (miRDB, miRTarBase, and TargetScan websites). Results Five IRMs were determined to develop the IRMRI. It can effectively segregate all HC cases from two different risk subgroups. We identified a marked discrepancy in survival outcome between the two groups by survival analysis and confirmed the reliability of IRMRI in two testing sets. Moreover, we collected 10 hub target genes (ESR1, IGF1, PDGFRB, JUN, MYC, ZWINT, MAD2L1, TOP2A, KIF11, and CDCA8) which were strongly linked to HC progression and malignant behavior. Conclusion We screened out five hub IRMs with clinical value and constructed a risk index model in HC, which can precisely assess the risk status and outcome of patients to a certain extent.


Introduction
Hepatocellular carcinoma (HC) is the sixth most common malignancy in the world, with more than 850,000 new cases each year [1]. HC ranks second among cancer-related causes of death and its incidence is increasing year by year. As a highly aggressive malignancy, HC kills 750,000 patients worldwide each year [2]. The incidence of HC has a large geographical heterogeneity, with approximately 85% of HC occurring in developing countries and regions and 72% of HC occurring in Asia [3]. Although there are many studies on the development of HC, it is still the tip of the iceberg for understanding its mechanisms. Currently, there are several well-defined risk factors for HC, including cirrhosis, virus infection, high alcohol intake, and aflatoxin B1. It is hard to diagnose HC given its insidious onset in early clinical stage [4]. Serum AFP test and imaging are the most common clinical tests, but these methods have limitations in the diagnosis of early stage. The survival outcome of HC is dismal owing to the tendency of metastasis and the unsatisfactory curative effect [5].
The immune microenvironment (IME) is a medium for the formation of immunocytes infiltrating in tumor tissues. Generated by tumor cells in their struggle with the immune system, IME is the condition and basis for tumor immune escape [6]. HC has unique self-protection mechanisms to evade host immune surveillance, such as secretion of immunosuppressive cytokines, abnormal expression of antigens, and alteration of the local IME [7]. For instance, TGF-β has a dual role in tumorigenesis. It blocks tumor cell viability and induces cell apoptosis in the early stage of disease, while it exerts an immunosuppressive role in the late stage of cancer. In HC, abnormal elevation of TGF-β1 suppresses the innate immune response and disrupts the antitumor immune response, which in turn facilitates tumor progression [8].
In our academic research, we unearthed the prognostic value of IRMs and set up a risk index of HC based on IRMs which could effectively forecast risk status and survival outcome of HC samples. In the future, personalized treatment for the HC will benefit from our constructed risk index.

Data Preparation.
In this project, the HC-associated dataset with RNA transcriptome data were collected from TCGA website. Based on the processing of the miRNA data, 375 HC cases were enrolled into our research. We also acquired 374 HC samples with complete mRNA data. In addition, the clinical traits and survival data were obtained from the UCSC Xena source. Each sample in the HC cohort with follow-up time < 30 days were removed. And a total of 245 IRMs were extracted from the Immune-miR database.

Determination of Differentially
Expressed miRNAs (DEMs) and mRNA (DEGs). Difference analysis was per-formed to collect DEMs and DEGs between HC and normal specimens by "limma" package in R [15]. We selected DEIRMs and DEGs on the basis of jlogFCj > 1 and adj-p < 0:05 as a filter.

Development of the Prognostic Index.
To set up an IRMrelated index (IRMRI), HC cases were equally and randomly divided into two sub-cohorts (training cohort and validation cohort). We first employed univariate analysis to screen possible IRMs showing dramatically prognostic ability in the training set. Then, the IRMs were processed by multivariate analysis to generate corresponding coefficients for IRMRI.
The risk value of IRMRI = ∑exp ðIRMsÞ * coef . The coef is the coefficient of each model IRMs.
2.4. Functional Enrichment Analysis. The target genes of miRNA were collected by three predictive online tools (miRDB, miRTarBase, and TargetScan websites). Next, the intersection target genes were collected by overlapping with the list of DEGs. R package "clusterProfiler" was instrumental in uncovering potential function of these intersection genes [16].

Identification of Hub
Gene. STRING is an online tool for examining protein-protein interactions (PPI) [17]. An interaction network was created by introduction of the target genes into STRING and visualized via Cytoscape tool [18]. Subsequently, we used CytoHubba algorithm to determine hub genes from the network based on degree score. At the same time, the prognostic power of target genes was detected by the Kaplan-Meier (KM) method.
2.6. Statistical Analysis. All statistical data was processed by R v.3.8.2. KM survival method contributed to assess the prognosis difference between two groups. The specificity and reliability of the index model were checked using ROC curves. Moreover, the Cox regressions were utilized to examine the independent value of clinical outcome.

Identification of Prognostic IRMs.
After processing prognosis data from the TCGA-HC, 184 HC cases were randomly assigned into the training set. In order to characterize the IRMs with prognostic value, univariate regression was conducted in the training set. We observed that a total of eight IRMs were tightly bound up with clinical outcome (Table 1).

Independence
Analysis of the IRMRI. We undertook an intensive analysis of clinical traits on age, gender, stage, and grading of HC cases to estimate the independent performance of the IRMRI. Univariate and multivariate methods both suggested that IRMRI was tightly related to patient outcome (Figures 4(a) and 4(b)). As indicated by Figure 4(c), the AUC of risk score was higher than other clinical factors, indicating a superior predictive power.
3.5. Functional Enrichment Analysis. There were 3,283 target genes determined by three miRNA-related databases. We performed further analysis by making predictions in three databases and taking target genes that were predicted in more than two databases. A total of 623 overlapping genes were then screened to detect the underlying function and pathway ( Figure 5(a)). In terms of GO, we observed that neurogenesis regulation, neuron differentiation, and urogenital development were greatly activated. As for KEGG, PI3K/ Akt signaling, MAPK signaling, and Rap1 signaling markedly enriched ( Figure 5(b)).

Discussion
HC is a classic inflammation-related cancer, and the IME plays a central part in the pathogenesis of HC [19]. IME is considered to be a key feature of cancer since the alterations in the IME are involved in all stages of malignant progression from the initial transformation stage to invasion and metastasis. Immunotherapy aims to provide more effective tumor cell targeting by enhancing existing tumor-specific immune responses [20].
In recent years, immunotherapy has been employed as an effective curative strategy for a variety of tumors, including HC [21]. In particular, therapies targeting immune checkpoints have achieved success and improved the clinical outcome of HC cases [22]. However, only a minority of patients benefit from immunotherapy due to the immunosuppressive status in IME [19]. Considering the prominence of IME in cancer progression, investigators should concentrate on uncovering new immune biomarkers and targets for HC management which can offer a reference for early diagnosis and prognosis determination.
Gene immunotherapy has become a promising approach for tumor treatment by restoring the function of tumor suppressor genes or stimulating the production of antitumor immune responses [23]. One of the effective medical methods is immunocytokine therapy which can be achieved by transfecting cytokines such as IL-2 directly in tumor and adjacent tissues [24]. Moreover, the clinical application of    Journal of Oncology immune checkpoint inhibitors opens up new mindsets for HC management. Immune checkpoints that have been extensively analyzed in relation to HC immune escape include PD-1/PD-L1 [25]. PD1 expressed on T cells, B cells, and natural killer (NK) cells could bind with PD-LI and PD-L2 ligands, inhibiting antigen-specific T cell activation and blocking the immune response of T cells in IME [22]. miRNAs have critical biological functions and their altered expression can contribute to cancer progression. Numerous reports have highlighted that miRNAs can regulate tumor initiation and progression as either pro-or anticancer factors [26]. The vast majority of HC can originate from cirrhosis of the liver due to various causes. As various etiologies lead to persistent liver injury and regeneration, individual HC etiologies also result in differential miRNA expression [27]. Hepatitis C virus infection is an integral factor in the pathogenesis of HC. Serum miRNA-27a may be used as an indicator of hepatitis C virus-induced HC [28]. In addition, Cao et al. found that upregulated miRNA-182-5p boosts HC initiation and progression [29]. Nevertheless, the regulatory role of IRMs in malignant behavior of HC needs a more exhaustive investigation.
In this work, an IRM-based prognostic index (has-miR-9-1, has-miR-30d, has-miR-139, has-miR-188, and has-miR-326) was created to analyze the risk status and survival outcome for HC cases. Our established IRMPI has turned out to display optimal independence with respect to survival outcome of patients. Survival analysis can dramatically differentiate the survival prognosis of the two groups. Meanwhile, we applied ROC analysis to detect the predictive reliability of IRMPI. Given the prominence of mRNA-miRNA interactions, we also determined the target genes of five model miR-NAs and obtained 10 hub genes (ESR1, IGF1, PDGFRB, Reviewing the previous literature, we found that these model miRNAs are more or less implicated in the formation of various tumors. As suggested by Wong CC et al., miR-139 could serve as a suppressor in HC and upregulation of miR-139 could inhibit the metastatic behavior of HC cell through suppression of ROCK2 [30]. NETA1 could increase the expression of TGF-β to facilitate HC cell viability through binding with miR-139 [31]. Also, miR-139 could mediate cell growth and metastasis in HC by targeting Wnt/TCF-4 axis [32]. As a promising noninvasive indicator in breast cancer (BC), miR-188 has been shown to regulate BC cell viability and metastasis by interacting with IL6ST [33]. In glioma, Li N et al. suggested that miR-188 might suppress cell cycle and cell development by inhibiting β-catenin [34].
Moreover, we found that among the 10 hub genes, 6 genes (CDCA8, ESR1, MAD2L1, KIF11, IGF1, and ZWINT) had significant prognostic value. Meanwhile, these 6 key genes are tightly bound up with tumor aggravation. For instance, Jeon et al. suggested that CDCA8 could regulate HC cell viability and stemness by targeting Akt/β-catenin signaling [35]. Aresti et al. detected the ERS1 expression in lung cancer and normal specimens and indicated its potential ability in evaluating survival outcome [36]. The function of MAD2L1 is related to cell mitosis. As revealed by Li et al., miR-200c-5p could block HC cell growth and survival by binding with MAD2L1 [37]. Similar to our results, Huang     [38]. Both KIF11 and ZWINT can affect HC progression by mediating cell proliferation [39,40]. The signaling pathway is a cellular reaction to signal reception and integration, which modulates gene expression and affects cell viability and apoptosis. Tumor-associated pathways engage in host of cellular activities and metabolic regulation. The activation of these pathways could contribute to the development of tumor malignant behavior and treatment failures [41]. Therefore, therapeutic strategies targeting signaling pathways have been the focus of cancer research, and various targeted inhibitors with promising benefits are becoming prevalent.
PI3K/Akt signaling pathway is mainly modulated by multiple genes such as PTEN, SHIP, and CTMP. PTEN blocks the dephosphorylation of PIP3 to PIP2, which could decrease the expression level of cellular PIP3 in cells and suppress the activation of Akt and its downstream molecules [42]. The downstream regulatory targets of the PI3K/Akt signaling include mTOR, Foxo, and GSK-3. Among them, the mTOR protein complex is a pivotal member. Akt can trigger mTOR kinase activity by negatively regulating the mTORC1 protein complex TSC1-TSC2, which mediates biological processes such as cell cycle, DNA damage repair, and glycogen synthesis [43]. In HC, Jiang et al. reported that PRMT9 could confer powerful migration ability to cancer cells by triggering Akt/GSK-3β pathway [44].
MAPK can be stimulated by mitogens, cytokines, and neurotransmitters to mediate cellular signals and exert biological effects by regulating cell growth, apoptosis, and autophagy. The MAPK pathway is composed of a conserved three-tier kinase pattern, including MKK, MKK, and MAPK, which can be triggered sequentially and act on downstream molecules, such as c-Jun and ATF2/6, to control the expression of specific genes and thus adjust the cell viability and differentiation [45].
Up to now, five parallel MAPK signaling pathways have been identified, including ERK1/2, SAPK, ERK5 p38, and JUK. Although each pathway is highly specific, in some cases, there is some crosstalk between them. ERK is mainly a transmitter of cell proliferation signals, while JNK and p38MAPK are mainly activated by various extracellular stimuli, causing a complex series of cellular stress transduction [46]. Several studies have revealed that MAPK signaling remains active at both transcriptional and translational levels in chemotherapy resistance HC cells [47]. Abnormal activation of MAPK pathway is responsible for the loss of cell differentiation and apoptosis, triggering abnormal cell proliferation and malignant transformation, and plays a central part in the mechanism of tumor drug resistance. In addition, upregulated ERK can promote cell division and proliferation by facilitating the cell cycle to enter S phase. Researchers have discovered that Ras mutations in cells induced abnormal upregulation of downstream MEK and ERK, resulting in abnormal proliferation of tumor cells [48].
In summary, we created a five IRM-based prognostic index which could play a central part in survival assessment for HC samples. In the future, our proposed index model might be a valuable informative tool offering clinical guide in HC management.

Data Availability
Public data were analyzed in this project. All data can be collected from TCGA database (https://portal.gdc.cancer.gov/).

Conflicts of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.