Profiling and Integrated Analysis of Differentially Expressed MicroRNAs as Novel Biomarkers of Hepatocellular Carcinoma

Hepatocellular carcinoma (HCC) is a heterogeneous disease that has multiple etiologies. It is the most common primary liver cancer, the sixth highest cause of cancer incidences, and the fourth highest cause of cancer-related deaths. The discovery of new biomarkers for the early detection, treatment, and prognosis of HCC would therefore be extremely useful. This study investigated differentially expressed ribonucleic acid (RNA) profiles by constructing a genome-wide profile of clinical samples. Differential expression analysis identified 1,280 differentially expressed messenger RNAs (dif-mRNAs), 99 differentially expressed microRNAs (dif-miRNAs), 181 differentially expressed long non-coding RNAs (dif-lncRNAs), and 31 differentially expressed circular RNAs (dif-circRNAs). Gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) path analysis were then conducted on these differentially expressed RNAs, revealing that they were clearly related to cell division, foreign body metabolism, and ribosome assembly. A competing endogenous RNA (ceRNA) network was then constructed based on the regulatory dif-miRNA-dif-mRNA and dif-miRNA-dif-lncRNA relationships. These results were also verified using HCC data from the Cancer Genome Atlas (TCGA); seven dif-miRNAs were verified in clinical samples by real-time quantitative polymerase chain reaction (RT-qPCR). Kaplan-Meier survival analysis revealed that the expression levels of Hsa-miR-1269a, Hsa-miR-421, and Hsa-miR-190b were correlated with overall survival. (P <0.05). Survival analysis of clinical samples showed that hsa-mir-1269a, hsa-mir-421 were associated with prognosis (p<0.05).This study revealed the general expression characteristics of specific differentially expressed miRNAs using a ceRNA network constructed from HCC samples. Hsa-mir-1269a, hsa-mir-421 may be promising candidates.

Introduction

Hepatocellular carcinoma (HCC) is the most common type of primary liver cancer; it is also one of the most common primary malignancies worldwide. According to global cancer statistics, HCC is also the fourth highest cause of cancer death (1–3). HCC usually develops in the context of chronic liver disease; its main causes are alcoholic and non-alcoholic fatty liver, and hepatitis B virus (HBV) and hepatitis C virus (HCV) infections (4). Therefore, it is urgent to better understand HCC tumor biology.

The rapid development of high-throughput next-generation sequencing (NGS) technologies has enabled the complete sequencing of entire genomes, permitting further analysis of the genomic profiles of some cancers (5–7). A recent ribonucleic acid (RNA) sequencing study of HCC found that a new type of vacuolar protein sorting 35 (VPS35) oncogene was significantly reduced in HCC cells. VPS35 plays a carcinogenic role by increasing the sorting and transport of fibroblast growth factor receptor 3 (FGFR3) (8). Recent studies have shown that competitive endogenous RNAs (ceRNAs) can act as natural decoys. The systematic functionalizations of non-coding RNAs (ncRNAs), pseudogenes, and circular RNAs (circRNAs) containing microRNA (miRNA) response elements (MREs) together form a complex miRNA-centered ceRNA network. This network constitutes a common microRNA library that affects gene expression and reveals new mechanisms of RNA interactions (9–12). As a ceRNA mechanism, the long-ncRNA (lncRNA) FAL1 can accelerate cell proliferation and metastasis by competitively binding to miR-1236 (13). However, few studies have used whole-transcriptome sequencing strategies to describe a transcriptional map, which could accurately detect the global gene expression profile.

Here, whole-transcriptome sequencing of tumors and corresponding adjacent non-tumors was conducted for six HCC patients. The expression levels of mRNA, lncRNA, miRNA, and circRNA were then analyzed in the cancer and paracancer groups, and their functional interactions were predicted. In addition, these results were validated using HCC data from the Cancer Genome Atlas (TCGA); a ceRNA network regulatory mechanism was then constructed. This study could offer a new molecular mechanism that could help to reveal the onset and progression of HCC, and its prognosis.

Materials and Methods

Patients and Tissue Collection

Cancerous and paracancerous tissue samples were obtained from six male patients with HCC, none of whom had received preoperative radiotherapy or chemotherapy. Patients were selected from the Department of Hepatobiliary and Pancreatic Surgery of the Affiliated Hospital of Qingdao University. All patients participating in the study signed an informed consent form. All specimens were examined for HCC by a pathologist prior to preservation (freezing via liquid nitrogen). This research was approved by the Medical Ethics Committee of the Affiliated Hospital of Qingdao University.

RNA Extraction and Quality Control

Total RNA was extracted from frozen tumor tissue and matched paracancerous tissue using TRIzol reagent (Thermofisher Science; Waltham, Massachusetts, USA). RNA degradation and contamination were monitored using 1.5% Sepharose gel, especially for deoxyribonucleic acid (DNA) contamination. RNA concentration and purity were measured using a NanoDrop 2000 spectrophotometer (ThermoFisher Scientific, Wilmington, DE). RNA integrity was assessed using an Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA) with the RNA Nano6000 assay kit.

Library Construction and Sequencing

RNA samples (2.5 ng) were used as input material to prepare libraries for RNA sequencing. Sequencing libraries were constructed and generated using the NEBNextR UltraTM small RNA Sample Prep Kit for Illumina R (NEB, USA), according to the manufacturer’s instructions. Index codes were added to attribute sequences to each individual sample. Finally, polymerase chain reaction (PCR) products were purified using the AMPure XP system (AMPure XP system); library quality was evaluated using the Agilent Bioanalyzer 2100 system. For clustering and sequencing, index-coded samples were clustered on the cBot cluster generation system, using the TruSeq PE Cluster Kitv3-cBot-HS (Illumina) according to the manufacturer’s instructions. After cluster generation, the library was sequenced on the Illumina platform and reads were generated.

Screening of Differently-Expressed (Dif)-Messenger RNAs (mRNAs), Dif-miRNAs, Dif-lncRNAs, and Dif-circRNAs

Samples with biological replicates were screened for differences using DEseq software; samples without biological replicates were screened for differences using EBseq software. Dif-mRNA, dif-circRNA, dif-lncRNA, and dif-miRNA were screened using fold changes greater than or equal to one and a false discovery rate (FDR) of <0.05 as screening criteria. Dif-mRNA, dif-lncRNA, and dif-miRNA were analyzed using the R package “clusterProfiler” for Gene Ontology-Biological Process (GO-BP) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment of differential genes (14–16).

Protein-Protein Interaction (PPI) Network and Module Analysis of Dif-mRNAs

STRING (version: 11.0; https://www.string-db.org/) was used to analyze the interaction relationships between differential proteins (17). The species was set to human and the import gene set was set to all degrees. To obtain the interaction pairs with the closest interaction relationship, the PPI arguments score was set to 0.99 (high confidence). The PPI network was constructed using Cytoscape software (version 3.8.0; https://cytoscape.org/) and the scores of each node were ranked to obtain the significant nodes. The most important cluster modules in the network were analyzed using Cytoscape’s plugin MCODE (18) (version 2.0.0; PPIMCODE https://apps.cytoscape.org/apps/). The screening threshold was set to 5.0 points, k-core=2, degree cutoff=2, node score cutoff=2, and max depth=100 for screening. GO-BP enrichment analysis was also performed for the important cluster module genes. A GO-BP threshold of P <0.05 was considered to imply a significant difference.

Construction of CeRNA Network

The ceRNA network was constructed using miRDB (V5.0; http://mirdb.org), miRTarBase (V7.0; https://maayanlab.cloud/Harmonizome/resource/MiRTarBase), miRWalk (http://mirwalk.umm.uni-heidelberg.de/) and TargetScan (V7.2; TSCAN.org/VERT_72/). These databases were searched for dif-miRNAs targeting dif-mRNAs to obtain dif-miRNA-dif-mRNA regulatory relationships (19). The dif-miRNA-dif-lncRNA regulatory relationships were obtained by searching for interactions between dif-lncRNAs and dif-miRNAs using the miRcode database (version 11; http://www.mircode.org/), RNAhybrid (https://bibiserv.cebitec.uni-bielefeld.de/rnahybrid/). The dif-miRNA-dif-mRNA and dif-miRNA-dif-lncRNA regulatory relationships were integrated and a composite lncRNA-miRNA-mRNA network was established using Cytoscape, according to the regulatory mechanism of the ceRNA network (20).

TCGA Liver Hepatocellular Carcinoma (LIHC) Database Validation

The data included 374 tumor specimens and 50 paracancerous normal specimens. The expression matrix and annotation information for mRNA, miRNA, and lncRNA were downloaded and the unpaired t-test in LIMMA software package was used to analyze the differential expression of tumor and normal tissues. A |logFC| value of >1 and a p value of <0.05 for mRNA and miRNA were used as cutoff points for the differential expression of lncRNAs. The differential expression results of TCGA data, and the dif-mRNAs, dif-miRNAs, and dif-lncRNAs obtained in this study were subjected to co-expression analysis, and a Venn graph was drawn.

Real-Time Quantitative PCR (RT-qPCR)

Thirty specimens of HCC tissues or normal liver tissues and histopathologically verified paracancerous tissues were selected from the Department of Hepatobiliary and Pancreatic Surgery, Affiliated Hospital of Qingdao University. Total RNA was extracted from HCC tissues and paraneoplastic tissues using Trizol (Takara). Total RNA was reverse transcribed into complementary DNA (cDNA) using a PrimeScript TM RT kit (Takara). RT-qPCR was performed on a LightCycler 480 using ChamQ Universal SYBR qPCR Master Mix (Vazyme) to determine miRNA expression (21). The study was approved by the ethics committee of Affiliated Hospital of Qingdao University.

Survival Analysis

Survival analysis of dif-miRNA using the Kaplan-Meier plotter (http://kmplot.com). The results were evaluated according to RT-qPCR. Grouping by 2^-ΔΔCT into high and low expression groups, survival analysis was performed according to hsa-miR-1269a, hsa-miR-421, and hsa-miR-190b.

Result

Dif-RNA Analysis

On the basis of the set screening conditions, a total of 1,280 dif-mRNAs were obtained, of which 904 were up-regulated and 376 were down-regulated. Furthermore, 99 dif-miRNAs were obtained, of which 31 were up-regulated and 68 were down-regulated; 181 dif-lncRNAs were obtained, of which 148 were up-regulated and 33 were down-regulated. Finally, 31 dif-circRNAs were obtained, of which 29 were up-regulated and two were down-regulated. The heat cluster maps of dif-mRNAs, dif-lncRNAs, dif-miRNAs, and dif-circRNAs showed that the tumor samples were clearly divided from the control samples (Figure 1), suggesting that the obtained differential expression analysis results were credible.

FIGURE 1

Figure 1 Heatmap of differential expression (N, Normal; T, Tumor) (A) differentially expressed mRNAs; (B) differentially expressed lncRNAs; (C) differentially expressed miRNAs; (D) differentially expressed circRNAs. Red indicates up-regulation and green indicates down-regulation.

Functional Enrichment Analysis of Dif-mRNAs

Following GO and KEGG pathway analysis, the top ten most enriched dif-mRNAs were identified. GO term analysis revealed that the genes displayed were related to biological processes such as cell division, and to the xenobiotic metabolism process (Figure 2A). KEGG analysis revealed that these genes were significantly related to fatty acid degradation, complement and coagulation cascades, and peroxisome signaling pathways, among other processes (Figure 2B).

FIGURE 2

Figure 2 Functional enrichment analysis of Dif-mRNAs: (A) Top 10 are enriched for Go-BP terms. (B) Top 10 enriched in KEGG pathway.

PPI Network and Module Extraction

The PPI network of dif-mRNAs consisted of 654 nodes and 2262 mutual pairs, among which nodes with high topological structure scores were considered to be the key nodes. Employing the Cytoscape plug-in MCODE, four sub-network modules (score R ≥5.0) were collected and extracted from the PPI network. Module A (score=26) contained 26 nodes and 325 interaction pairs, module B (score = 10.333) contained 19 nodes and 93 interaction pairs, module C (score = 6.8) contained 16 nodes and 51 interaction pairs, and module D (score=5.273) contained 12 nodes and 29 reciprocity pairs (Figure 3A). In addition, Go-BP enrichment analysis of the interacting genes in these modules revealed the top five enrichment results of each module; these were selected for display in order of importance. Among them, the genes of module A were found to be closely related to organelle fission and nuclear division, and the genes of module B were found to be related to platelet degranulation and humoral immune response. The genes in module C were revealed to be related to the negative regulation of the cell cycle process, and to the cell cycle G1/S phase transition. The genes in module D, meanwhile, were found to be related to the G1/S phase transitions of the cell cycle and the mitotic cell cycle (Figure 3B).

FIGURE 3

Figure 3 Four modules extracted from the PPI network and GO-BP enrichment analysis of the four modules (A) The four modules in the ppi network with a score>5; (B) Top5 of GO-BP enrichment analysis of genes in the four modules.

Enrichment Analysis of Dif-lncRNA, Dif-miRNA Targeted Regulation Gene

Functional enrichment analysis of dif-lncRNAs (divided into Cis and Trans dif-lncRNAs) target gene and dif-miRNAs target gene was performed, with the top ten richest results being selected. In GO-BP analysis, Cis dif-lncRNA target gene was found to be significantly related to protein k63-linked ubiquitination and the activation of NF-Kappa B-inducing Kinase activity, among other processes. (Figure 4A). Tran dif-lncRNA target gene was shown to be associated with processes such as cell division and the viral modulation of host morphology or physiology (Figure 4B). Dif-miRNA target gene was revealed to be significantly related to axon guidance and the neurotrophic tyrosine kinase (Trk) receptor signaling pathway, among other processes (Figure 4C). KEGG analysis revealed that Cis dif-lncRNA target gene was significantly correlated with retinol metabolism and chemical carcinogenesis (Figure 4D) whereas Tran dif-lncRNA target gene was significantly correlated with the cell cycle and RNA transport (Figure 4E). Dif-miRNA target gene was revealed to be associated with pathways in cancer, aldosterone synthesis, and secretion (Figure 4F).

FIGURE 4

Figure 4 Dif-lncRNAs、Dif-miRNAs GO-BP and KEGG enrichment analysis (A) Top10 GO-BP enrichment analysis of cis-lncRNAs target genes; (B) Top10 GO-BP enrichment analysis of tran-lncRNAs target genes; (C)Top10 GO-BP enrichment analysis of miRNAs target genes; (D) Top10 KEGG pathway enrichment analysis of cis-lncRNAs target genes; (E) Top10 KEGG pathway enrichment analysis of tran-lncRNAs target genes; (F) Top10 KEGG pathway enrichment analysis of miRNAs target genes.

CeRNA Network Construction

According to the mutual regulation formation and information of the dif-miRNA-dif-mRNA and dif-miRNA-dif-lncRNA relationships, lncRNAs and mRNAs that were remarkable differentially expressed and mutually regulated on alike miRNAs were screened. Cytoscape was used to construct the lncRNA-miRNA-mRNA ceRNA network (Figure 5).

FIGURE 5

Figure 5 Construction of ceRNA networks by differentially expressed mRNAs, miRNAs and lncRNAs.

TCGA LIHC Data Verification

After setting the screening thresholds at p <0.05 and fold change >1, screening the TCGA LIHC dif-mRNA, dif-lncRNA, and dif-miRNA data revealed 4,840 dif-mRNAs, 2,616 dif-lncRNAs, and 251 dif-miRNAs. These were compared with the dif-mRNAs, dif-lncRNAs, and dif-miRNAs obtained from Seq analysis and the dif-lncRNAs and dif-miRNAs identified during co-expression analysis. Thus, 805 dif-mRNAs, 6 dif-lncRNAs, and 7 dif-miRNAs were screened; the resulting vine is shown in Figures 6A–C. Among these RNAs, all seven dif-miRNAs (Table 1) were up-regulated in HCC tissues.

FIGURE 6

Figure 6 TCGA HCC database and Seq differentially expressed RNAs co-expression analysis (A) Dif-mRNAs co-expression analysis; (B) Dif-lncRNAs co-expression; (C) Dif-miRNAs co-expression.

TABLE 1

Table 1 Sequence of each miRNA.

RT-qPCR and Survival Analysis

Seven dif-miRNAs were detected by qRT-PCR, all of which showed significantly higher expression levels in all HCC tissues (Figures 7A–G), consistent with the results obtained from the preliminary data. Survival analysis of these seven dif-miRNAs in the Kaplan-Meier database (http://kmplot.com) showed that the overall survival rates of HCC patients with high expression of has-miR-1269a, has-miR-421, and has-miR-190b were significantly lower than those of patients with lower expression of these targets (p<0.05 Figures 7H–N). Patients were divided into two groups according to RT-qPCR 2^-ΔΔCT: low expression group and high expression group. Subsequent Kaplan-Meier analysis showed that high expression of has-miR-1269a, has-miR-421 was associated with shorter overall survival (p<0.05, Figures 8A–C).

FIGURE 7

Figure 7 Dif-miRNAs expression in tissue and K-M survival analysis of overall survival (OS) rate (A) has-miR-1269a Up-regulation in HCC (*P < 0.05); (B) has-miR-421 Up-regulation in HCC (**P < 0.01); (C) has-miR-4326 Up-regulation in HCC (*P < 0.05); (D) has-miR-7706 Up-regulation in HCC (*P<0.05); (E) has-miR-944 Up-regulation in HCC (**P < 0.01); (F) has-miR-190b Up-regulation in HCC (P < 0.05); (G) has-miR-217 Up-regulation in HCC (*P < 0.05); (H) Correlation of hsa-miR-1269a with OS (p < 0.05); (I) Correlation of hsa-miR-421 with OS (*p < 0.05); (J) Correlation of hsa-miR-4326 with OS; (K) Correlation of hsa-miR-7706 with OS; (L) Correlation of hsa-miR-944 with OS; (M) Correlation of hsa-miR-190b with OS (p < 0.05); (N) Correlation of hsa-miR-217 with OS.

FIGURE 8

Figure 8 K-M survival analysis of overall survival (OS) of hsa-miR-1269a, hsa-miR-421, and hsa-miR-190b in clinical samples. (A) Correlation of hsa-miR-1269a with OS (p < 0.05) (B) Correlation of hsa-miR-421 with OS (p < 0.05) (C) Correlation of hsa-miR-4326 with OS.

Discussion

HCC is one of the most common malignant tumors of the gastrointestinal tract; it is the sixth-most common malignant tumor worldwide, and has the fourth highest mortality rate. Although surgical resection, radiotherapy, chemotherapy, and liver transplantation have all improved significantly in recent years (22, 23), the prognosis remains unsatisfactory. The recurrence rate remains high five years after surgery (24). Therefore, elucidating the molecular mechanisms and processes of HCC is important for discovering new therapeutic targets and improving the clinical prognoses of patients (13).

CeRNA is a usually used to refer to a class of RNAs, incorporating lncRNAs, miRNAs, mRNAs, and circRNAs. These substances can act as sponges after competitively binding to shared miRNAs. CeRNA networks link protein-coding mRNAs to ncRNAs (e.g., miRNAs and lncRNAs) after participating in the development of tumorigenesis (25, 26).

Here, whole transcriptome sequencing of cancer and paraneoplastic tissues from six HCC patients resulted in the identification of 1,280 dif-mRNAs, 99 dif-miRNAs, 181 dif-lncRNAs and 31 dif-circRNAs. Functional enrichment analysis of dif-mRNAs showed that these differential genes were majorly involved in cell division, xenobiotic metabolic processes, and ribosome assembly. A PPI network was constructed, hub genes were analyzed and identified, and GO-BP enrichment analysis was performed for four modules. These four modules were significantly associated with organelle fission, platelet degranulation, the negative regulation of cell cycle processes, and cell cycle G1/S phase transition. GO-BP and KEGG enrichment analyses were performed for dif-miRNAs, dif-lncRNAs target gene. A ceRNA network of lncRNA-miRNA-mRNA was constructed based on the mutual regulatory dif-miRNA-dif-mRNA and dif-miRNA-dif-lncRNA relationships. In addition, most dif-mRNAs and dif-lncRNAs, and a few dif-miRNAs, were successfully verified by TCGA data. The differences highlighted in both sequencing and TCGA data might be related to sample differences or threshold selection. RT-qPCR validation and survival analysis (using the Kaplan-Meier database) was performed for the seven screened dif-miRNAs.

Few studies have used ceRNA networks to assess the prognosis of HCC, and the discovery of potential biomarkers is crucial to improve the diagnoses and prognoses of HCC patients. It has been shown that the ceRNA networks of breast cancer tissues and matched normal tissues differ, with some ceRNAs being activated in cancerous tissues, leading to the development of breast cancer; they are inactive in normal cells, however (27). Those ceRNAs that are differently expressed in cancer and normal tissues could be used as potential markers. It is well known that miRNAs are an integral part of the cancer development process (28).Thus, studies into the prognostic relevance of miRNA regulatory aspects are essential. In this context, here seven co-expressed dif-miRNAs were validated in HCC tissues, consistent with transcriptome sequencing results and differential miRNA expression in TCGA HCC. Then, survival analysis was performed, revealing that Hsa-miR-1269a, Hsa-miR-421, and Hsa-miR-190b were significantly correlated with overall survival (p <0.05). It has been found that Hsa-miR-1269a is upregulated in advanced colorectal cancer, and that it forms a regenerative feedback loop with the transforming growth factor-beta (TGF-β) signaling pathway to promote metastasis in colorectal cancer (29). Hsa-miR-1269a has been shown to be highly expressed in uveal melanoma and clear cell renal cell carcinoma, interrelated with poor patient prognosis. Therefore, it can be used as a prognostic effect to predict patient survival (30, 31). In serum exosomes, Hsa-miR-1269a plays an oncogenic role in non-small cell lung cancer (NSCLC); it can be used as a diagnostic marker (32). Hsa-miR-421 has been revealed to be highly expressed in NSCLC. Furthermore, upon binding to kelch-like ECH-associated protein 1 (KEAP1) three prime untranslated region (3’UTR), it has been used to forecast low survival in NSCLC (33). On the contrary, the expression of Hsa-miR-421 expression has been found to be downregulated in breast cancer tissues and metastatic cell lines. Decreased levels of its expression have been observed to be associated with lymph node metastasis, recurrence, metastasis, or TNM staging (34). The low expression of Hsa-miR-190b in breast cancer, meanwhile, has been associated with better prognoses (35). Although dif-miRNAs have been validated to some extent, the interconnections of dif-miRNAs in prediction networks have not been discussed. As previously reported, miRNA-1269a facilitates the proliferation and apoptosis of glioma cells by directly targeting ATRX (36). The lncRNA DLEUI increases the expression of rho associated coiled-coil containing protein kinase 1 (ROCK1) via Hsa-miR-421, thereby promoting the progress of papillary thyroid carcinoma (37). Furthermore, the lncRNA tumor suppressor candidate 8 (TUSC8) can inhibit the evolution of papillary thyroid carcinoma via the miR-190b myosin regulatory light chain interacting protein (MYLIP) axis, which in turn can inhibit breast cancer growth and metastasis (38). Here, miRNA-1269a was found to target binding to six transmembrane epithelial antigen of prostate 4 (STEAP4) and kinase insert domain receptor (KDR). STEAP4 is associated with adipocyte metabolism and mediates hepatocellular carcinogenesis, whereas KDR is a major growth factor of endothelial cells. It is the main mediator in inducing endothelial cell proliferation, survival, migration and tubule morphogenesis; it is also closely related to tumorigenesis and metastasis.

This study also has some limitations, however. First, it is limited by its small sample size regarding the supply of the general situation of the HCC transcriptome; this may have led to analytical bias. In addition, the functions of the identified dif-miRNAs were not been investigated in depth. The predicted lncRNA-miRNA-mRNA interaction networks need to be further confirmed by performing in vitro cellular function experiments, and by using animal models to understand the functions and mechanisms behind them. Researchers are currently investigating the direct functions of identified dif-miRNAs and the molecular mechanisms of targeting RNAs.

Conclusions

In this study, a ceRNA network was constructed based on whole transcriptome sequencing data. When screened in combination with TCGA, seven differentially expressed miRNAs were identified; further validation revealed that three miRNAs Hsa-miR-1269a, Hsa-miR-421, were significantly associated with prognosis. Therefore, these differential miRNAs are expected to be potential biomarkers or therapeutic targets for the prognosis of HCC.

Data Availability Statement

The data presented in the study are deposited in the CNGB Sequence Archive (CNSA) of China National GeneBank DataBase (CNGBdb) repository, accession number CNP0002592.

Ethics Statement

The studies involving human participants were reviewed and approved by Medical Ethics Committee of the Affiliated Hospital of Qingdao University. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

YX, YW, and WX performed the experiments, designed the research, analyzed the data, and wrote the manuscript. HZ, KunL, KuiL, and WZ analyzed the data, and wrote the manuscript. CZ and JC provided reagents and intellectual guidance. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the Taishan Scholars Program of Shandong Province (grant number 2019010668), and the Shandong Higher Education Young Science and Technology Support Program (grant number 2020KJL005).

Conflict 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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Yang JD, Hainaut P, Gores GJ, Amadou A, Plymoth A, Roberts LR. A Global View of Hepatocellular Carcinoma: Trends, Risk, Prevention and Management. Nature Reviews. Gastroenterol Hepatol (2019) 16(10):589–604. doi: 10.1038/s41575-019-0186-y

CrossRef Full Text | Google Scholar

2. Global Burden of Disease Cancer Collaboration, Fitzmaurice C, Allen C, Barber RM, Barregard L, Bhutta ZA, et al. Global, Regional, and National Cancer Incidence, Mortality, Years of Life Lost, Years Lived With Disability, and Disability-Adjusted Life-Years for 32 Cancer Groups, 1990 to 2015: A Systematic Analysis for the Global Burden of Disease Study. JAMA Oncol (2017) 3(4):524–48. doi: 10.1001/jamaoncol.2016.5688

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Siegel RL, Miller KD, Jemal A. Cancer Statistics, 2020. CA: Cancer J Clin (2020) 70(1):7–30. doi: 10.3322/caac.21590

PubMed Abstract | CrossRef Full Text | Google Scholar

4. El-Serag HB, Rudolph KL. Hepatocellular Carcinoma: Epidemiology and Molecular Carcinogenesis. Gastroenterology (2007) 132(7):2557–76. doi: 10.1053/j.gastro.2007.04.061

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Ma YS, Huang T, Zhong XM, Zhang HW, Cong XL, Xu H, et al. Proteogenomic Characterization and Comprehensive Integrative Genomic Analysis of Human Colorectal Cancer Liver Metastasis. Mol Cancer (2018) 17(1):139. doi: 10.1186/s12943-018-0890-1

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Ma YS, Wu ZJ, Zhang HW, Cai B, Huang T, Long HD, et al. Dual Regulatory Mechanisms of Expression and Mutation Involving Metabolism-Related Genes FDFT1 and UQCR5 During CLM. Mol Ther Oncolytics (2019) 14:172–8. doi: 10.1016/j.omto.2019.04.008

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Harding JJ, Nandakumar S, Armenia J, Khalil DN, Albano M, Ly M, et al. Prospective Genotyping of Hepatocellular Carcinoma: Clinical Implications of Next-Generation Sequencing for Matching Patients to Targeted and Immune Therapies. Clin Cancer Res (2019) 25(7):2116–26. doi: 10.1158/1078-0432.CCR-18-2293

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Zhang G, Tang X, Liang L, Zhang W, Li D, Li X, et al. DNA and RNA Sequencing Identified a Novel Oncogene VPS35 in Liver Hepatocellular Carcinoma. Oncogene (2020) 39(16):3229–44. doi: 10.1038/s41388-020-1215-6

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Yu F, Liu JB, Wu ZJ, Xie WT, Zhong XJ, Hou LK, et al. Tumor Suppressive microRNA-124a Inhibits Stemness and Enhances Gefitinib Sensitivity of Non-Small Cell Lung Cancer Cells by Targeting Ubiquitin-Specific Protease 14. Cancer Lett (2018) 427:74–84. doi: 10.1016/j.canlet.2018.04.022

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Ma YS, Lv ZW, Yu F, Chang ZY, Cong XL, Zhong XM, et al. MicroRNA-302a/D Inhibits the Self-Renewal Capability and Cell Cycle Entry of Liver Cancer Stem Cells by Targeting the E2F7/AKT Axis. J Exp Clin Cancer Res: CR (2018) 37(1):252. doi: 10.1186/s13046-018-0927-8

CrossRef Full Text | Google Scholar

11. Sumazin P. Et Al, An Extensive microRNA-Mediated Network of RNA-RNA Interactions Regulates Established Oncogenic Pathways in Glioblastoma. Cell (2011) 147(2):370–81. doi: 10.1016/j.cell.2011.09.041

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Tay Y, Kats L, Salmena L, Weiss D, Tan SM, Ala U, et al. Coding-Independent Regulation of the Tumor Suppressor PTEN by Competing Endogenous mRNAs. Cell (2011) 147(2):344–57. doi: 10.1016/j.cell.2011.09.029

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Long J, Bai Y, Yang X, Lin J, Yang X, Wang D, et al. Construction and Comprehensive Analysis of a ceRNA Network to Reveal Potential Prognostic Biomarkers for Hepatocellular Carcinoma. Cancer Cell Int (2019) 19:90. doi: 10.1186/s12935-019-0817-y

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Ju Q, Zhao YJ, Ma S, Li XM, Zhang H, Zhang SQ, et al. Genome-Wide Analysis of Prognostic-Related lncRNAs, miRNAs and mRNAs Forming a Competing Endogenous RNA Network in Lung Squamous Cell Carcinoma. J Cancer Res Clin Oncol (2020) 146(7):1711–23. doi: 10.1007/s00432-020-03224-8

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Zhao Y, Zhang H, Ju Q, Li X, Zheng Y. Comprehensive Analysis of Survival-Related lncRNAs, miRNAs, and mRNAs Forming a Competing Endogenous RNA Network in Gastric Cancer. Front Genet (2021) 12:610501. doi: 10.3389/fgene.2021.610501

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Ju Q, Zhao YJ, Dong Y, Cheng C, Zhang S, Yang Y, et al. Identification of a miRNA-mRNA Network Associated With Lymph Node Metastasis in Colorectal Cancer. Oncol Lett (2019) 18(2):1179–88. doi: 10.3892/ol.2019.10460

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Zhang YJ, Ma YS, Xia Q, Yu F, Lv ZW, Jia CY, et al. MicroRNA−mRNA Integrated Analysis Based on a Case of Well−Differentiated Thyroid Cancer With Both Metastasis and Metastatic Recurrence. Oncol Rep (2018) 40(6):3803–11. doi: 10.3892/or.2018.6739

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Bandettini WP, Kellman P, Mancini C, Booker OJ, Vasu S, Leung SW, et al. MultiContrast Delayed Enhancement (MCODE) Improves Detection of Subendocardial Myocardial Infarction by Late Gadolinium Enhancement Cardiovascular Magnetic Resonance: A Clinical Validation Study. J Cardiovasc Magn Reson (2012) 14:83. doi: 10.1186/1532-429X-14-83

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Du JW, Li GQ, Li YS, Qiu XG. Identification of Prognostic Biomarkers Related to the Tumor Microenvironment in Thyroid Carcinoma. Sci Rep (2021) 11(1):16239. doi: 10.1038/s41598-021-90538-3

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Smoot ME, Ono K, Ruscheinski J, Wang PL, Ideker T. Cytoscape 2.8: New Features for Data Integration and Network Visualization. Bioinformatics (2011) 27(3):431–2. doi: 10.1093/bioinformatics/btq675

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Xu X, Kong X, Liu T, Zhou L, Wu J, Fu J, et al. Metastasis-Associated Protein 1, Modulated by miR-30c, Promotes Endometrial Cancer Progression Through AKT/mTOR/4E-BP1 Pathway. Gynecol Oncol (2019) 154(1):207–17. doi: 10.1016/j.ygyno.2019.04.005

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Tang A, Hallouch O, Chernyak V, Kamaya A, Sirlin CB. Epidemiology of Hepatocellular Carcinoma: Target Population for Surveillance and Diagnosis. Abdominal Radiol (2018) 43(1):13–25. doi: 10.1007/s00261-017-1209-1

CrossRef Full Text | Google Scholar

23. Rudnick SR, Russo MW. Liver Transplantation Beyond or Downstaging Within the Milan Criteria for Hepatocellular Carcinoma. Expert Rev Gastroenterol Hepatol (2018) 12(3):265–75. doi: 10.1080/17474124.2018.1417035

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Li J, Wu PW, Zhou Y, Dai B, Zhang PF, Zhang YH, et al. Rage Induces Hepatocellular Carcinoma Proliferation and Sorafenib Resistance by Modulating Autophagy. Cell Death Dis (2018) 9(2):225. doi: 10.1038/s41419-018-0329-z

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA Hypothesis: The Rosetta Stone of a Hidden RNA Language? Cell (2011) 146(3):353–8. doi: 10.1016/j.cell.2011.07.014

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Qi X, Zhang DH, Wu N, Xiao JH, Wang X, Ma W. ceRNA in Cancer: Possible Functions and Clinical Implications. J Med Genet (2015) 52(10):710–8. doi: 10.1136/jmedgenet-2015-103334

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Conte F, Fiscon G, Chiara M, Colombo T, Farina L, Paci P. Role of the Long Non-Coding RNA PVT1 in the Dysregulation of the ceRNA-ceRNA Network in Human Breast Cancer. PLoS One (2017) 12(2):e0171661. doi: 10.1371/journal.pone.0171661

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Patel N, Garikapati KR, Makani VKK, Nair AD, Vangara N, Bhadra U, et al. Regulating BMI1 Expression via miRNAs Promote Mesenchymal to Epithelial Transition (MET) and Sensitizes Breast Cancer Cell to Chemotherapeutic Drug. PLoS One (2018) 13(2):e0190245. doi: 10.1371/journal.pone.0190245

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Bu P, Wang L, Chen KY, Rakhilin N, Sun J, Closa A, et al. miR-1269 Promotes Metastasis and Forms a Positive Feedback Loop With TGF-β. Nat Commun (2015) 6:6879. doi: 10.1038/ncomms7879

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Zhao D-D, Zhao X, Li W-T. Identification of Differentially Expressed Metastatic Genes and Their Signatures to Predict the Overall Survival of Uveal Melanoma Patients by Bioinformatics Analysis. Int J Ophthalmol (2020) 13(7):1046–53. doi: 10.18240/ijo.2020.07.05

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Qin S, Shi X, Wang C, Jin P, Ma F. Transcription Factor and miRNA Interplays Can Manifest the Survival of ccRCC Patients. Cancers (2019) 11(11). doi: 10.3390/cancers11111668

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Wang X, Jiang X, Li J, Wang J, Binang H, Shi S, et al. Serum Exosomal miR-1269a Serves as a Diagnostic Marker and Plays an Oncogenic Role in Non-Small Cell Lung Cancer. Thorac Cancer (2020) 11(12):3436–47. doi: 10.1111/1759-7714.13644

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Duan FG, Wang MF, Cao YB, Dan Li, Li RZ, Fan XX, et al. MicroRNA-421 Confers Paclitaxel Resistance by Binding to the KEAP1 3'UTR and Predicts Poor Survival in Non-Small Cell Lung Cancer. Cell Death Dis (2019) 10(11):821. doi: 10.1038/s41419-019-2031-1

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Pan Y, Jiao G, Wang C, Yang J, Yang W. MicroRNA-421 Inhibits Breast Cancer Metastasis by Targeting Metastasis Associated 1. Biomed pharmacother = Biomed Pharmacother (2016) 83:1398–406. doi: 10.1016/j.biopha.2016.08.058

CrossRef Full Text | Google Scholar

35. Gong PJ, Shao YC, Huang SR, Zeng YF, Yuan XN, Xu JJ, et al. Hypoxia-Associated Prognostic Markers and Competing Endogenous RNA Co-Expression Networks in Breast Cancer. Front Oncol (2020) 10:579868. doi: 10.3389/fonc.2020.579868

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Zhang Y, Wang Q, Luo N, Liu J, Ren H, Shao X, et al. MicroRNA-1269a Promotes Proliferation and Arrest of Apoptosis of Glioma Cells by Directly Targeting ATRX. Front Oncol (2020) 10:563901. doi: 10.3389/fonc.2020.563901

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Li R, Wan T, Qu J, Yu Y, Zheng R. Long Non-Coding RNA DLEUI Promotes Papillary Thyroid Carcinoma Progression by Sponging miR-421 and Increasing ROCK1 Expression. Aging (2020) 12(20):20127–38. doi: 10.18632/aging.103642

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Zhao L, Zhou Y, Zhao Y, Li Q, Zhou J, Mao Y. Long Non-Coding RNA TUSC8 Inhibits Breast Cancer Growth and Metastasis via miR-190b-5p/MYLIP Axis. Aging (2020) 12(3):2974–91. doi: 10.18632/aging.102791

PubMed Abstract | CrossRef Full Text | Google Scholar