Comprehensive Analysis of Competitive Endogenous Rnas Network Associated With Hepatocellular Carcinoma

Background: Increasing evidences show that long non-coding RNA (lncRNA) plays the role of competitive endogenous RNAs (ceRNAs) in the development and progression of cancers. The purpose of our study was to identify potential lncRNA biomarkers that serve as a therapeutic target and prognostic biomarker in HCC. Methods: The differential expression of RNAs was examined using the edgeR package. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were used to predict potential functions. Survival analysis and ROC curve analysis were performed to predict ceRNA network had signicant prognostic value. The biological functions of MYLK-AS1 on HCC cells were studied by RNA interference approaches in vitro. Cell proliferation, migration and invasion were detected by Cell Counting Kit-8(CCK-8) assay, wound healing and transwell assay relatively. The mechanism of competitive endogenous RNAs (ceRNAs) were predicted and veried by bioinformatic analysis, western blot analysis and luciferase assays. Results: The newly constructed ceRNA network comprised 76 HCC specic lncRNAs, 15 miRNAs, and 35 mRNAs from the database (Targetscan, miRTarBase, and miRDB). 10 differentially expressed lncRNAs and 10 differentially expressed mRNAs were signicantly associated with overall survival in HCC (P value < 0.05). Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes(KEGG) pathways enrichment analysis results showed the differentially expressed mRNAs were involved primarily in the most signicant cancer-related signaling pathways. ROC curve analysis demonstrated that MYLK-AS1— has-mir-424—CCNE1 ceRNA network had signicant prognostic value. The expression of MYLK-AS1 was up-regulated in HCC cells and tissues. Biological function analyses indicated that down-regulation of MYLK-AS1 suppressed cell proliferation, migration and invasion . MYLK-AS1 and miR-424-5p bound directly and reversibly to each other. MYLK-AS1 could positively regulate CCNE1 expression by competitively binding to miR-424-5p. Conclusion: study into the lncRNA-related ceRNA network in HCC and the MYLK-AS1 may be a candidate biomarker for molecular diagnosis and prognosis monitoring of HCC.


Introduction
Hepatocellular carcinoma (HCC) is a kind of malignant tumor with high morbidity and mortality [1][2] . More than 70% of new cases of HCC occur in Asia each year, and more than 50% of these new cases occur in China [3] . The treatment and adverse reactions of HCC have always been global challenges. Therefore, the search for early diagnostic markers and more effective and safer treatments are signi cance for improving the clinical strategies and prognosis of HCC. LncRNAs are transcripts of more than 200 nucleotides that have been shown to be involved in a variety of biological processes such as chromosomal silencing [4] , chromatin modi cation [5] , transcriptional activation [6] and transcriptional interference [7] . There is continuous evidence that lncRNAs can participate in a variety of physiological and pathological processes in human cancer, affecting cell proliferation, migration and invasion. For example, LncRNA-HOST2 can promote cell proliferation, migration and invasion and inhibit cell apoptosis in human HCC cell line SMMC-7721 [8] ; LncRNA WWOX-AS1 inhibits the proliferation, migration and invasion of osteosarcoma cells [9] ; Lnc-MMP2-2 might regulate the migration and invasion of lung cancer cells into the vasculature by promoting MMP2 expression, suggesting this lncRNA as a novel therapeutic target and predictive marker of tumor metastasis in lung cancer [10] ; Long non-coding RNA lnc-GNAT1-1 inhibits gastric cancer cell proliferation and invasion through the Wnt/β-catenin pathway in Helicobacter pylori infection [11] .
LncRNA acts as a competitive endogenous RNA (ceRNA) that binds to microRNAs (miRNAs) to regulate its downstream target genes [12][13][14][15] . At the same time, miRNA is an important regulator of HCC tumorigenesis and development. miR-199a-3p regulates MTOR and PAK4 pathways and inhibits tumor growth in hepatocellular carcinoma, miR-199a-3p may be promising as an HCC treatment option [16] . APA Zekri, A. N. et al. identi ed a miRNA panel comprised of four miRNAs (miR-192, miR-122, miR-181b and miR-125a-5p) that may serve as a molecular tool for characterization of the CD133+ cells associated with different stages of hepatocarinogensis [17] . MiR-18a may serve as a prognostic biomarker of HCC as it is demonstrated to carry out a decisive role in HCC progression by promoting HCC cell invasion, migration, and proliferation through targeting Bcl2L10 [18] .
In present research was based on TCGA, we analyzed the RNA expression pro les from 374 tumors and 50 adjacent normal tissues in HCC. As a result, 1082 differentially expressed lncRNAs (DElncRNAs), 122 differentially expressed miRNAs (DEmiRNAs) and 1992 differentially expressed mRNAs (DEmRNAs) were identi ed in our study. In order to elucidate the interactions and valid potential crosstalk between RNAs, we successfully established the HCC associated ceRNA network based on bioinformatics generated from Targetscan, miRTarBase, and miRDB, which included 76 lncRNAs, 15 miRNAs and 35mRNAs.
Furthermore, genes of ceRNA network were analyzed for overall survival to identify prognostic genes with clinical characteristics. In addition, we identi ed MYLK-AS1 regulated tumorigenesis of HCC. We found that MYLK-AS1 was markedly upregulated in HCC cells. Knocking down MYLK-AS1 could inhibit HCC cell proliferation and migration via functioning as a ceRNA for miR-424-5p, thereby preventing its association with CCNE1. Collectively, the results showed that MYLK-AS1 is a novel tumor biomarker that can be used as a potential target for clinical diagnosis and treatment of HCC.

Materials And Methods
Patients and samples from the TCGA database To search the speci c of lncRNAs, miRNAs and mRNAs in HCC, we downloaded gene expression data from the TCGA database (https://portal.gdc.cancer.gov/, up to April 9, 2020, including 374 HCC tissues and 50 adjacent normal tissues).
RNA sequence data processing and differential expression analysis The differentially expressed mRNAs (DEmRNAs), lncRNAs (DElncRNAs) and miRNAs (DEmiRNAs) between the HCC tissues and adjacent normal tissues was identi ed by EdgeR package in R (version 3.4.1), and P < 0.01 and |logFC| > 2 were set as the cut-off criteria. The ggplot2 packages in R were used to visualize volcano plots.
Functional enrichment analysis GO enrichment analysis were performed with the Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/), P<0.05 was set as the cut-off criterion. KEGG pathway analysis were performed with packages in R, the threshold was P<0.05, R clusterPro ler package were used to predict potential functions for the DEmRNAs in the ceRNA network.

Human HCCtissues and adjacent normal tissues microarray
The expression of MYLK-AS1 was detected by human HCC tissue and adjacent normal tissue microarrays (Lnc cDNA-HLivH090Su01, OUTDO BIOTECH, Shanghai, China). These tissue microarrays included 64 cases and 26 HCC tissues and adjacent normal tissues, respectively. The Taizhou Hospital Ethics Committee of Zhejiang Province authorized all experiments in patient organizations in this study and obtained informed consent from each patient participating in the study.

Cell culture
Human liver cancer cell lines (Huh7 and HepG2) and normal liver epithelial cell line (LO2) were obtained from the Chinese Cell Bank of the Chinese Academy of Sciences (Shanghai, China). All cells were cultivated in Dulbecco modi ed Eagle medium (DMEM) together with 10% fetal bovine serum (FBS; Gibco Waltham, MA), 100U/mL penicillin as well as 100μg/mL streptomycin (both from Sigma-Aldrich, St Louis, MO) in a humidi ed temperature at 37°C with 5% CO2 RNA extraction and real-time PCR analysis Total RNAs were extracted from cells using RNAiso Plus (TaKaRa, Tokyo, Japan). A reverse transcription kit (RR036A; TaKaRa) was used to transcribe total RNA and produce complementary DNA (cDNA). Realtime PCR was performed with SYBR Premix Ex Taq II (TaKaRa) and NanoDrop 2000c (Thermo Scienti c, Waltham, MA, USA). Expression of genes was normalized to that of glyceraldehyde 3-phosphate dehydrogenase (GAPDH) or U6. Relative RNA expressions were measured with ABI Prism 7500 Software v2.0.6 (Thermo Fisher Scienti c, Waltham, MA) and calculated based on the 2 −ΔΔC t method.

Western blot analysis
Proteins were isolated from tissues by lysing frozen tissues in radioimmunoprecipitation assay (RIPA) buffer (Sigma-Aldrich). Equal amounts of protein (50 μg) were separated by 12% sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS PAGE) and transferred to a polyvinylidene uoride (PVDF) membrane (Thermo Fisher Scienti c, Billerica, MA). The membrane was blocked with 5% bovine serum albumin (BSA) overnight and incubated with primary antibodies including CCNE1 and GAPDH (Cell Signaling Technology). GAPDH was used as an internal control. Protein bands were visualized by the enhanced chemiluminescence (ECL) detection system (Applygen Technologies, Beijing, China).

Cell proliferation assay
The cell suspension (100 μL/well) was seeded in a 96-well plates. The plate was pre-incubated in an incubator (at 37 °C, 5% CO 2 ). Add 10 μL of CCK-8 solution to each well. Incubate the plates in the incubator for 24, 48, 72, 96 hours. The absorbance at 450 nm was measured with a microplate reader.

Colony formation assay
Huh7 and HepG2 cells were trypsinized to make a single cell suspension, 500 cells were inoculated into each well and the medium was changed every three days. Culture was stopped after 3 weeks, the cells were washed with phosphate-buffered saline (PBS), xed with methanol and stained with crystal violet (Sinopharm Chemical Reagent, Beijing, China). The colonies were counted under a uorescence microscope (Olympus Corporation, Tokyo, Japan).

Wound healing assay
The marker pen draws a horizontal line evenly behind the 6-well plates (Corning), with a horizontal line every 0.5 ~ 1.0 cm. 5 × 10 5 cells were seeded into to the well. Cell scratch line perpendicular to the horizontal line was made using a 200 μL pipette tip. 6-well plates were incubated in the incubator(at 37°C , 5% CO 2 ). The sample was taken out and photographed at 0, 24 hours.
Transwell invasion assay 5 × 10 4 cells/200 ul of serum-free cell suspension was added to the Matrigel Invasion Chambers (size 8 μm; BD Biosciences, Franklin Lakes, NJ, USA), 600 μL of medium containing 10% FBS was added to the lower chamber of 24-well plates (Corning). Continue to cultivate for 24 hours, the chamber was removed and the phosphate-buffered saline (PBS) were used to wash the cells twice, the cells were xed in methanol, and the cells were stained with crystal violet. The cells were observed, and were counted 5 elds of view under the microscope (Olympus, Tokyo, Japan).

Dual luciferase reporter assay
Wild plasmids was constructed by cloning into the psi-CHECK-2 vector (Promega, Madison, WI, USA) using XhoI and NotI sites. Mutant plasmids were generated using the QuickChange® Site-Directed Mutagenesis kit (Stratagene, Agilent Technologies, Wilmington, DE, USA). The appropriate plasmid and miR-424-5p mimic or mimic control were co-transfected into HEK293T cells (1.0 × 10 5 ), and the luciferase assay was evaluated 48 hours after transfection using the Dual-luciferase Reporter Assay System (Promega). Renilla luciferase activity was used for standardization.

Statistical analysis
The clinical data on the patients were combined with HCC data in TCGA to evaluate the prognostic value of differential RNAs in the ceRNA network. Survival curves were generated using the survival package in R for samples with differentially expressed mRNAs and lncRNAs. All results presented as the mean ± standard deviation (S.D.). Oneway analysis of variance (ANOVA) was used to examine statistical comparisons among groups. Two-tailed Student's t-test was performed to compare the means of values from different experiments. Differences were considered statistically signi cant if P < 0.05. All statistical analyses were performed using GraphPad Prism 5.0 (GraphPadSoftware, Inc., La Jolla, CA, USA). Each experiment was performed for three times.

Predictions of mRNAs and lncRNAs targeted by miRNAs
Next, we predicted the mRNAs and lncRNAs that were targeted by miRNAs, focusing on the relationship the 122 differentially expressed miRNAs and 1082 differentially expressed lncRNAs above. Only 15 of 122 differentially expressed miRNAs were predicted to target 76 of 1082 differentially expressed lncRNAs (table1). The relationships between these 15 differentially expressed lncRNA-targeting miRNAs were used to predict the targeted mRNAs using Targetscan, miRTarBase, and miRDB. Then, 15 HCC-speci c miRNAs were predicted to target the 35 mRNAs (table2).  To reveal how lncRNA may mediate transcription in HCC by affecting mRNA and miRNA binding, a ceRNA network based on the lists of DElncRNAs, DEmiRNAs, and DEmRNAs was constructed and visualized using Cytoscape software. As shown in Fig. 2, the lncRNA-miRNA-mRNA network was comprised of 15 miRNA nodes, 35 mRNA nodes, 76lncRNA nodes.

Functional enrichment analysis based on mRNAs
To establish context of ceRNA network, we inferred the roles of each lncRNA based on the functions of connected mRNAs. lncRNAs were typically central and connected to one or more mRNAs in the network.
The top 16 highly enriched GO terms of biological process (BP), cellular component (CC) and molecular function (MF). GO analysis was conducted to ascertain the signaling cascade that the 25 genes participate (Fig. 5). Finally, KEGG pathway analysis revealed that 9 pathways were signi cantly enriched, particularly mRNAs involved in cancer (Fig. 6).

Differential expression of MYLK-AS1 in HCC cells and tissues
In our study, Our results showed that MYLK-AS1 expression was signi cantly increased in two HCC cell lines (Huh7, HepG2) compared with that in normal HCC epithelial cell line LO2 (Fig.8A).MYLK-AS1 expression was higher in HCC tissues than adjacent normal tissues by human HCC tissue and adjacent normal tissue microarrays. (Fig.8B, C). In addition, MYLK-AS1 expression was further analysed according to patients' clinical pathological features. Higher expression of MYLK-AS1 correlated with larger tumour size, advanced TNM stage (Fig.8D, E). The area under the ROC curve of MYLK-AS1 was 0.75 (Fig.8F). Kaplan-Meier survival analysis and log-rank tests showed that higher MYLK-AS1 expression was associated with shorter survival time (Fig. 8G). Collectively, these ndings indicate that MYLK-AS1 is a potential biomarker for diagnosis and prognosis in HCC.

MYLK-AS1 silencing inhibits HCC cells proliferation,migration and invasion
To examine the biological functions of MYLK-AS1, the si-MYLK-AS1 1#, si-MYLK-AS1 2#, si-MYLK-AS1 3#, or si negative control (NC) were separately transfected into Huh7 and HepG2 cells. We found that MYLK-AS1 expression was remarkably reduced in si-MYLK-AS1 3# transfected into HCC cells compared with si NC (Fig. 9A). The CCK8 assay and plate clone formation assay was used to determine the role of MYLK-AS1 in cell growth. In Huh7 and HepG2 cells, MYLK-AS1 knockdown led to reduced cell proliferation compared with that observed in the si NC cells (Fig. 9B, C). Furthermore, in the wound healing assay, knockdown of MYLK-AS1 contributed to slower scratch healing (Fig. 9D). Transwell assay revealed that the invasion ability of cells in which MYLK-AS1 was silenced was suppressed compared with that of si NC (Fig.9E). These results reported that MYLK-AS1 could promote the proliferation and metastasis of HCC cells.

MYLK-AS1 functions as a ceRNA and sponges miR-424-5p in HCC cells.
To investigate the molecular mechanism by which MYLK-AS1 acting as ceRNAs. First of all, MYLK-AS1 knockdown signi cantly increased the expression level of miR-424-5p (Fig. 10A), This nding suggested that MYLK-AS1 may function as a ceRNA of miRNAs. To determine this hypothesis, We then performed dual luciferase reporter assays to con rm the prediction analysis. HEK293T cells were transfected with a luciferase plasmid harboring the sequence of MYLK-AS1 together with plasmids encoding the miRNAs or a control sequence. We found that miR-424-5p could suppress MYLK-AS1-driven luciferase activity, and the suppression ability of miR-424-5p is stronger (Fig. 10B). To determine whether miR-424-5p functions as a tumor suppressor in HCC cells, we transfected Huh7 and HepG2 cells with miR-424-5p mimic or inhibitor (Fig. 10C). Real-time PCR analysis con rmed that the expression of MYLK-AS1 was lower in miR-424-5p mimic compared with mimic control, the expression of MYLK-AS1 was higher in miR-424-5p inhibitor compared with inhibitor control (Fig. 10D). We then performed CCK-8 and colony formation found that cell proliferation and colony formation ability were signi cantly reduced by overexpression of miR-424-5p (Fig. 10E,F)and signi cantly enhanced by silencing of miR-424-5p expression (Fig.11A,B). the wound healing assay and transwell assay found that cell migration and invasive ability were signi cantly reduced by overexpression of miR-424-5p (Fig. 10G,H) and signi cantly enhanced by silencing of miR-424-5p expression (Fig.11C,D) CCNE1 is a miR-424-5p target gene and is indirectly regulated by MYLK-AS1 To determine the ceRNA network between MYLK-AS1, miR-424-5p and its targets in HCC. We found that knockdown of MYLK-AS1 also signi cantly reduced CCNE1 mRNA and protein levels in Huh7 and HepG2 cells (Fig. 12A). To determine whether CCNE1 is regulated by miR-424-5p liver cancer cells, we measured CCNE1 mRNA and protein levels when miR-424-5p was over-expressed or inhibited in Huh7 and HepG2 cells. We found that CCNE1 mRNA and protein levels were signi cantly decreased by miR-424-5p overexpression (Fig. 12B), in contrast, CCNE1 mRNA and protein levels were signi cantly increased by miR-424-5p inhibition (Fig. 12C). Nest, we performed luciferase reporter assays driven by the wild-type 3' UTR sequence of CCNE1, which contains the predicted miR-424-5p binding site (wt-CCNE1), or mutant constructs containing a mutation in the miR-424-5p-binding sites (mut-CCNE1). These plasmids were cotransfected into HEK293T cells together with miRNA mimic control or miR-424-5p mimic. The results showed that wt-CCNE1-driven luciferase expression was signi cantly reduced by co-transfection with the miR-424-5p mimic compared with the control, but this repression was abolished by mutation of the putative miR-424-5p-binding site in the CCNE1 3'UTR (Fig. 12D). Taken together, these results indicate that miR-424-5p regulates CCNE1 expression in liver cancer cells by directly binding to the predicted site in the 3' UTR of CCNE1 mRNA.
miR-424-5p plays a role in the relationship between MYLK-AS1and CCNE1 To determine whether miR-424-5p is involved in mediating the effects of MYLK-AS1 in HCC cells, Huh7 and HepG2 cells were co-transfected with si-MYLK-AS1 3# and miR-424-5p inhibitor. Notably, In addition, proliferation and colony forming assays revealed that inhibition of the miR-424-5p promoted the proliferation of Huh7 and HepG2 cells, and this effect was partly reversed by co-transfection with si-MYLK-AS1 3# (Fig. 13A and B), the wound healing assay and transwell invasion assay found that was partially rescued by co-transfection with miR-424-5p inhibitor( Fig. 13C and D). To determine whether miR-424-5p plays a role in the relationship between si-MYLK-AS1 3# and CCNE1, we examined cells cotransfected with si-MYLK-AS1 3# and the miR-424-5p inhibitor. Indeed, the suppression of CCNE1 protein levels induced by si-MYLK-AS1 3# was effectively reversed by the miR-424-5p inhibitor (Fig. 13E).
Collectively, these data suggest that MYLK-AS1 modulates the expression of CCNE1 by posttranscriptional regulation of miR-424-5p.

Discussion
A large amount of evidence has demonstrated that ceRNA plays an important role of human diseases, including thyroid cancer [20] , non-small cell lung cancer [21] , and ovarian cancer [22] . However, the role and mechanism of ceRNA in HCC remains unclear.
In this study, RNA sequencing data and the clinical characteristics of HCC were obtained from the Cancer Genome Atlas database. DElncRNAs, DEmRNAs, and DEmiRNAs were identi ed between HCC and normal HCC tissue samples. Subsequently, the lncRNA-miRNA-mRNA ceRNA network of HCC was established. We identi ed a novel lncRNA-MYLK-AS1 associated with HCC, which is signi cantly upregulated in HCC cell lines. This research demonstrated that MYLK-AS1 modulate cell proliferation, migration and invasion in HCC cells. Luciferase reporter assay con rmed that miR-424-5p was a target of MYLK-AS1 in HCC cells. These ndings indicated that MYLK-AS1 played an oncogenic role in HCC and could be considered as a potential prognostic indicator for HCC.
Although lncRNA has received extensive attention in recent years, miRNAs also deserves more attention. Undoubtedly, the cancer-related signaling pathways based on miRNAs regulation is indispensable.
Misregulated expression miRNAs is reported to play various roles in tumorigenesis. It is reported that miR-424-5p functions as tumor suppressor gene in many cancers. For example, the inhibition of has-mir-424 in SNHG12-depleted cells partially reversed the effects on cervical cancer cell apoptosis, adhesion and invasion [23] . Another report shows that has-mir-424 -SMAD7 pathway contributed to ESCC invasion and metastasis and up-regulation of has-mir-424 perhaps provided a strategy for preventing tumor invasion, metastasis [24] . In this study, we found that miR-424-5p was signi cantly down-regulated in HCC cells, and miR-424-5p mimic or inhibitor could impair or promote HCC cell proliferation, migration and invasion. Our ndings uncover the signi cance of the interaction between MYLK-AS1 and miR-424-5p in tumorigenesis given that MYLK-AS1 exerts oncogenic behavior partly via sponging miR-424-5p in HCC cells. In addition, miR-424-5p functions as a tumor suppressor and directly targets CCNE1 as potential prognostic markers in epithelial ovarian cancer [25] .
The miRNAs-targeted genes involved in ceRNA network were performed the GO term analysis and KEGG pathways enrichment analysis. These genes, E2F2, E2F7, E2F1,CCNE1, CDC25A and KIF23 not only enriched in GO term analysis or KEGG pathway but also included in the ceRNA network. We found that these genes play a key regulatory role in the occurrence and development of tumors [26][27][28][29][30] . Few studies have reported the role of CCNE1 in HCC, so the molecular mechanism of CCNE1 in HCC was also our focus. We conducted luciferase reporter assays and veri ed that miR-424-5p targeted CCNE1 mRNA at its 3 UTR. Moreover, miR-424-5p mimic inhibited CCNE1 protein expression, and miR-424-5p inhibitor promoted CCNE1 protein expression. Simultaneously, knockdown of MYLK-AS1 also signi cantly reduced CCNE1 protein level. Thus, we con rmed CCNE1 as the direct target of miR-424-5p. MYLK-AS1 may act as a ceRNA to regulate the expression of CCNE1 by repressing its inhibitor miR-424-5p.

Conclusion
In summary, we found that some lncRNAs and mRNAs were signi cantly associated with overall survival in HCC patients. Importantly, we successfully constructed the lncRNA -related ceRNA network, bringing new approach to lncRNA research in HCC, and providing novel lncRNA MYLK-AS1 as candidate prognosis biomarkers or potential therapeutic targets.   The most signi cant ten lncRNAs associated with overall survival in ceRNA network in HCC on Kaplan-Meier survival curves. The horizontal axis represents overall survival time (years), vertical axis represents survival function. (P < 0.05).

Figure 4
The most signi cant ten protein-coding genes associated with overall survival in ceRNA network in HCC on Kaplan-Meier survival curves. The horizontal axis represents overall survival time (years), vertical axis represents survival function. (P < 0.05).