Construction of a Myc-associated ceRNA network reveals a prognostic signature in hepatocellular carcinoma.

Hepatocellular carcinoma (HCC) remains an extremely lethal disease worldwide. High-throughput methods have revealed global transcriptome dysregulation; however, a comprehensive investigation of the complexity and behavioral characteristics of the competing endogenous RNA (ceRNA) network in HCC is lacking. In this study, we extracted the transcriptome (RNA) sequencing data of 371 HCC patients from The Cancer Genome Atlas platform. With the comparison of the high Myc expression (Mychigh) tumor and low Myc expression (Myclow) tumor groups in HCC, we identified 1,125 differentially expressed (DE) mRNAs, 589 long non-coding RNAs (lncRNAs), and 93 microRNAs (miRNAs). DE RNAs predicted the interactions necessary to construct an associated Myc ceRNA network, including 19 DE lncRNAs, 5 miRNAs, and 72 mRNAs. We identified a significant signature (long intergenic non-protein-coding [LINC] RNA 2691 [LINC02691] and LINC02499) that effectively predicted overall survival and had protective effects. The target genes of microRNA (miR)-212-3p predicted to intersect with DE mRNAs included SEC14-like protein 2 (SEC14L2) and solute carrier family 6 member 1 (SLC6A1), which were strongly correlated with survival and prognosis. With the use of the lncRNA-miRNA-mRNA axis, we constructed a ceRNA network containing four lncRNAs (LINC02691, LINC02499, LINC01354, and NAV2 antisense RNA 4), one miRNA (miR-212-3p), and two mRNAs (SEC14L2 and SLC6A1). Overall, we successfully constructed a mutually regulated ceRNA network and identified potential precision-targeted therapies and prognostic biomarkers.

Introduction

Liver cancer is categorized as primary hepatic cancer and metastatic hepatic cancer.1, 2, 3 Hepatocellular carcinoma (HCC) is the fifth-most common cancer worldwide and is a highly lethal malignancy, accounting for nearly 80% of all primary liver malignancies.4, 5, 6 According to 2018 Global Cancer statistics, there have been more than 100,000 new cases and deaths worldwide every year.

Although there have been great advancements in cancer therapy, the 5-year overall survival (OS) rate for HCC is still only 12.10%.8, 9, 10, 11 Surgical resection, liver transplantation, and chemotherapy are common therapies,12, 13, 14, 15 but they are appropriate only for patients with early-stage disease., Therefore, it is essential to understand the mechanisms underlying the pathogenesis of HCC and identify novel biomarkers.

In the past several years, the advancement of high-throughput technologies has provided more opportunities for biomarker identification.18, 19, 20, 21 Non-coding RNAs (ncRNAs), which include long ncRNAs (lncRNAs) and microRNAs (miRNAs), play a role in cancer progression and are potential biomarkers and therapeutic targets.22, 23, 24, 25 lncRNAs are transcripts with more than 200 nucleotides, which previously were thought to have no biological function. But in fact, lncRNAs are involved in several biological processes, including cell differentiation, cancer proliferation, and metastasis.26, 27, 28, 29 Competing endogenous RNA (ceRNA), including lncRNA, can interact with mRNA by competitively combining different miRNAs.30, 31, 32 miRNAs are a class of ncRNA molecules of 21–24 nucleotides that mediate negative post-transcriptional regulation.33, 34, 35 When the miRNA arm-imbalance mechanism is broken in the cell, dysregulation of downstream tumor-suppressor genes or oncogenes controlled by aberrant miRNAs occurs, leading to cancer development.36, 37, 38, 39, 40 For example, lncRNA LPP antisense RNA 2 has carcinogenic effects and promotes cell proliferation and metastasis through microRNA (miR)-7-5. The well-studied tumorigenic lncRNA, which is highly upregulated in liver cancer (HULC), serves as a ceRNA network through miR-372. AGAP2 antisense RNA 1, a competitive lncRNA, functions as an oncogene and upregulates annexin A11 expression through miR-16-5p, promoting proliferative capacity in liver cancer. These ceRNA networks provide a novel perspective and offer insight into undetected biomarkers for the early diagnosis and treatment of cancer.

The Myc gene is located on chromosome 8q24.21, which is dysregulated in most human neoplasia. This region is frequently genetically amplified in various human cancers. Moreover, MYC is overexpressed in more than 50% of tumors, including HCC. In this study, we first used the median expression level of MYC to divide the 371 HCC samples into two groups for subsequent analysis: high Myc expression (Mychigh) tumor and low Myc expression (Myclow) tumor. We used the transcriptome sequencing data of 371 lncRNAs, 371 mRNAs, and 367 miRNAs of HCC tumor tissue and adjacent normal tissue from The Cancer Genome Atlas (TCGA) platform to identify key differentially expressed (DE) RNAs, and we constructed an associated Myc ceRNA network to reveal the underlying mechanism in HCC carcinogenesis. We identified a prognostic signature (long intergenic non-protein-coding [LINC] RNA 2691 [LINC02691] and LINC02499) that effectively predicts OS and has protective effects.

We constructed a critical ceRNA network that included four lncRNAs (LINC02691, LINC02499, LINC01354, and NAV2 antisense RNA 4 [NAV2-AS4]), one miRNA (miR-212-3p), and two mRNAs (SEC14-like protein 2 [SEC14L2] and solute carrier family 6 member 1 [SLC6A1]). This study provided a better understanding of the pathogenesis of liver cancer from the perspective of polygenic association, thus offering novel insights into targeted combination therapies.

Results

Study process of transcriptome data

The study flow diagram is shown in Figure 1. We divided the transcriptome data of 371 HCC tissues into Mychigh tumor and Myclow tumor groups based on the standard median expression of Myc. We selected DE lncRNAs, miRNAs, and mRNAs from the Mychigh tumor group and Myclow tumor group. We considered lncRNA-miRNA to be a true interaction target by miRcode. We explored the miRDB database and TargetScan for miRNA-mRNA target prediction to construct a ceRNA network, followed by analyses of expression and survival.

Figure 1

Research diagram of the Myc ceRNA network in HCC

Clinical analysis of Myc overexpression in HCC

To explore whether Myc affects gene expression in liver cancer, we divided liver cancer patients into Mychigh tumor and Myclow tumor groups based on the median values of Myc in this study. The mRNA and protein levels of Myc were presented in various normal organs in the Human Protein Atlas (HPA) database (Figures 2A and S1A). The expression of Myc was high in cancer cell lines (Figure S1B). The region of Myc that is genomically altered in liver cancer was expressed through amplification and was found to be altered in 66 (18%) of 366 patients in HCC (Figures 2B and 2C). Myc expression was higher in HCC tissues (n = 369) compared with in normal liver tissues (n = 160, p < 0.01; Figure 2D). The expression of Myc increased as tumor invasion worsened (Figure 2E). Kaplan-Meier survival analysis showed the prognostic potential of Myc expression. These data suggested that Myc was upregulated in liver cancer samples.

Figure 2

Functional analysis of Myc

(A) Transcription levels of Myc in various normal organs using the HPA database. (B and C) Frequency of genetic alterations associated with Myc. (D) Expression of Myc was analyzed between cancer tissues (n = 369) and normal tissues (n = 160) in HCC. (E) Distribution of Myc expression between TNM stages.

DE RNAs in HCC

To better understand the relationship between DE RNAs and tumorigenesis-associated HCC, we downloaded gene-expression microarrays from TCGA database to search for DE RNAs. The gene-expression microarrays contained 371 samples (lncRNAs, mRNAs); 367 samples (miRNAs) were defined as Mychigh and Myclow based on the median expression levels of Myc. We analyzed the DE genes between the Mychigh and Myclow tumor groups with the standard of p < 0.05 and fold change (FC) ≥ 1.5. In total, we found 589 lncRNAs, 93 miRNAs, and 1,125 mRNAs of DE RNAs among the groups. The upregulated DE RNAs included 222 (38%) lncRNAs, 65 (70%) miRNAs, and 696 (62%) mRNAs. The downregulated DE RNAs included 367 (62%) lncRNAs, 38 (30%) miRNAs, and 429 (38%) mRNAs. In Figure 3, the volcano plots show the distribution of DE RNAs (Figure 3A), and the heatmap describes 15 significant DE RNAs (Figure 3B).

Figure 3

Identification of DE genes

(A) Volcano maps of DE lncRNAs, miRNAs, and mRNAs between two groups: Mychigh and Myclow in HCC. (B) Heatmap of DE lncRNAs (top), miRNAs (middle), and mRNAs (bottom). Red and blue spots represent significantly upregulated and downregulated RNAs, respectively.

Functional enrichment analysis of DE mRNAs

To better comprehend the mechanisms involved in HCC, we explored the function of 1,125 DE mRNAs from Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses using the Metascape database (Figures 4A and 4B). The most enriched GO terms in cellular component (CC), biological process (BP), and molecular function (MF) were “side of membrane,” “response to toxic substance,” and “organic anion transmembrane transporter activity,” respectively. We found the oxidation-reduction process, peroxisome proliferator-activated receptor (PPAR) signaling, and peroxisome pathway to participate in HCC, according to KEGG pathway enrichment analysis.

Figure 4

Functional enrichment analysis of DE mRNAs

(A) KEGG pathway of DE mRNAs. (B) MF, BP, and CC of DE mRNAs.

Construction of the ceRNA network and identification of hub RNAs

We used miRNA target informatic tools to identify DE RNAs in the ceRNA. The ceRNA network was constructed using Cytoscape 3.7 (Figure 5A) and included 19 lncRNAs, 5 miRNAs, and 72 mRNAs. The ceRNA network included 19 lncRNAs (11 downregulated and 8 upregulated) and 5 miRNAs (4 downregulated and 1 upregulated). In addition, 72 mRNAs (58 upregulated and 14 downregulated) and 5 miRNAs (4 downregulated and 1 upregulated) were used to construct the ceRNA network. We employed a large ceRNA network based on node-weighting arithmetic to recognize highly interacted hub clustering (Figure 5B). The hub RNAs contained 14 lncRNAs (6 upregulated and 8 downregulated), 5 miRNAs (2 upregulated and 3 downregulated), and 11 mRNAs (9 upregulated and 2 downregulated).

Figure 5

Myc ceRNA network in HCC

(A) Blue rhombus represents downregulated lncRNAs, red ellipses represent upregulated miRNAs, and blue triangles represent downregulated mRNAs; red rhombus represents upregulated lncRNAs, blue ellipses represent downregulated miRNAs, and red triangles represent upregulated mRNAs. (B) The hub 30 genes in the ceRNA network.

Comprehensive analysis of DE lncRNAs in the ceRNA network

lncRNAs have an impact on the responses of miRNAs and mRNAs, commanding the upstream portion of the ceRNA network. The expression of lncRNAs is associated with OS in cancer patients.46, 47, 48, 49 Therefore, we analyzed 19 DE lncRNAs of ceRNA in this study and found differential expression in the Mychigh tumor group (n = 185) and Myclow tumor group (n = 186) using the heatmap (Figure 6A). Among the 19 DE lncRNAs, we screened potential prognosis-associated lncRNAs by univariate and multivariate Cox regression analyses. The final risk score method established that an lncRNA signal was a significant prognostic factor for HCC, as follows: prognostic index (PI) = (−0.6621 × expression level of LINC02691) + (−0.2047 × expression level of LINC02499). Notably, the classification of these two groups using the risk score improved the distinguishing ability and predictive power regarding OS based on Kaplan-Meier (p = 0.0055) and time-dependent receiver operating characteristic (ROC) curve analyses (area under the curve [AUC] = 0.625). LINC02691 and LINC02499 were found to be protective factors for HCC (Figure 6B). lncRNAs of the ceRNA network play a role in the cytoplasm and are involved in the regulation of downstream genes. Hence, we used the lncLocator online tool to predict the location of LINC02691, LINC02499, LINC01354, and NAV2-AS4. We found that the subcellular locations of LINC02691 did not possess a sequences’ record in LNCipedia, which was not predicted. The transcripts analysis of LINC02499, LINC01354, and NAV2-AS4 obtained from LNCipedia indicated their location in the cytoplasm using lncLocator online prediction (Figure 6C).

Figure 6

Comprehensive analysis of DE lncRNAs in the ceRNA network

(A) Differential expression of 19 lncRNAs in the ceRNA network by a heatmap. (B) Validation of the two lncRNA signatures by Kaplan-Meier survival and ROC curves based on risk scores. (C) Prediction of lncRNA subcellular localization; a score ranging from 0 and 1 was given. (D) Expression of LINC02691, LINC02499, LINC01354, and NAV2-AS4 was analyzed in the Mychigh tumor group (n = 185) and Myclow tumor group (n = 186) in HCC. (E) Survival analysis.

To identify potentially significant RNAs regarding the survival of patients with HCC, we conducted DE analysis, clinical-stage analysis, and Kaplan-Meier survival curve for each RNA in the ceRNA network. A comparison of consequences follows: 4 lncRNAs (LINC02691, LINC02499, LINC01354, and NAV2-AS4) were suppressively expressed in the Mychigh tumor group (n = 185) relative to the Myclow tumor group (n = 186; Figure 6D); 4 lncRNAs were strongly correlated with OS (p < 0.05) (Figure 6E); and the remaining 15 lncRNAs were not significantly correlated with OS in HCC (Figure S2A).

We next explored the relationship of four lncRNAs according to clinical pathological stage. Following worsening tumor invasion, the expression of the four lncRNAs decreased (Figure S2B). As further verification of this finding, the difference in LINC02691, LINC02499, LINC01354, and NAV2-AS4 expression in liver cancer tissues (n = 371) was statistically significant (p < 0.05) compared with normal liver tissues (n = 50) (Figure S2C). The 50 paired cancer and paracancer tissues were employed to validate the results of DE lncRNAs (Figure S2D).

Functional analysis of miR-212-3p as a target

We searched for potentially significant miRNA in HCC. Five DE miRNAs (4 downregulated and 1 upregulated [miR-212-3p]) in the ceRNA network were associated with OS. We constructed a Kaplan-Meier survival curve for four miRNAs (Figure 7A). Moreover, patients with high miR-212-3p expression had a poor prognosis compared with those with low miR-212-3p expression. The results showed that miR-212-3p was a factor that promoted tumor growth. miR-212-3p was highly expressed in the Mychigh tumor group (n = 184) compared with the Myclow tumor group (n = 18; Figure 7B). The expression of miR-212-3p increased with high tumor-node-metastasis (TNM) stage in 371 liver cancers (Figure 7C).

Figure 7

Expression and survival analysis of DE miRNAs

(A) Kaplan-Meier survival analysis of miR-212-3p, miR-217, miR-216b-5p, miR-375, and miR-146b-5p in 371 HCC tissues. (B) The miR-212-3p expressed in the Mychigh tumor group (n = 184) and Myclow tumor group (n = 183). (C) The miR-212-3p expression in the TNM stage in 371 liver cancer samples. (D) Expression level of miR-212-3p in 371 HCC tissues compared with 50 adjacent normal tissues. (E) In 50 pairs of HCC tissues, miR-212-3p was highly expressed in liver cancer. (F) A schematic diagram of the miR-212-3p of 3′ UTR in four lncRNA-binding sites. (G) The high-risk group had high miRNA expression and low lncRNA expression. The low-risk group had low miRNA expression and high lncRNA expression.

We performed further validation with cancer and paracancer data and compared adjacent normal tissues. We found that the expression level of miR-212-3p was significantly upregulated in 371 HCC tissues. In 50 paired HCC tissues, miR-212-3p was highly expressed in liver cancer relative to nontumor tissues (Figures 7D and 7E). In addition, through the informatics websites, we found that miR-212-3p included specific linking of the 3′ UTR region of the lncRNA gene. The 3′ UTR binding location of miR-212-3p of four lncRNAs is shown in Figure 7F. We combined miRNA and lncRNA to verify the ceRNA hypothesis. The high-risk group had high miRNA expression and low lncRNA expression, whereas the low-risk group had low miRNA expression and high lncRNA expression. We conducted survival analysis comparing the high-risk group with the low-risk group (Figure 7G).

Comprehensive analysis of DE mRNAs in the ceRNA network

To further understand the biological functions, we explored the prognostic values of the 72 DE mRNAs involved in the ceRNA network. In the ceRNA network, we found that SEC14L2 and SLC6A1 were downstream target genes of miR-212-3p. According to Kaplan-Meier survival analysis, SEC14L2 and SLC6A1 were significantly associated with prolonged prognosis (p < 0.05), whereas the remaining the mRNAs in the ceRNA network had no survival significance in HCC (Figure S3). SEC14L2 and SLC6A1 were overexpressed in the Myclow tumor group and were associated with low-grade TNM stage (Figures 8A and 8B). Further validation showed that the expression of SEC14L2 and SLC6A in 371 liver cancers was lower than that in 50 paracancer samples; 50 samples were paired to explain the expression of the genes in cancer and paracancer tissues (Figures 8C and 8D). No genetic alterations were found in SEC14L2 and SLC6A1 in HCC; however, SEC14L2 was altered in 1 (0.3%) of 366 patients, and SLC6A1 was altered in 5 (1.4%) of 366 patients with liver cancer in the cBioPortal database. In addition, SEC14L2 and SLC6A1 had a lower expression in liver cancer (n = 10) compared with paired-normal tissue samples (n = 10) taken from the Gene Expression Omnibus (GEO) profiles (Figures 8E and 8F).

Figure 8

Analysis of SEC14L2 and SLC6A1

(A and B) Comparison of SEC14L2 and SLC6A1 expression in neoplasm development and association with prolonged prognosis. (C and D) Verification of SEC14L2 and SLC6A1 expression in normal tissues and HCC tissues in liver. (E and F) Verification of SEC14L2 and SLC6A1 in HCC relative to paired normal tissue samples from GEO profiles. (G) Immune cell infiltration of SEC14L2 and SLC6A in HCC using the TIMER database.

Expression of SEC14L2 and SLC6A1 in various cancers

We analyzed the transcription and protein levels of SEC14L2 and SLC6A1 in various organ tissues using the HPA database (Figures S4 and S5). To further evaluate SEC14L2 and SLC6A1 expression in human cancer, we used RNA sequencing data to examine SEC14L2 and SLC6A1 expression between tumor and normal tissues. The expression of SEC14L2 was significantly lower in some cancers, including bladder urothelial carcinoma, cholangiocarcinoma, kidney chromophobe, lung adenocarcinoma, and HCC. SLC6A1 expression was lower in most cancers, except renal cell clear carcinoma and HCC (Figure S6). The expression of SEC14L2 and SLC6A1 was low in various cancer cell lines (Figure S7).

Association of immune infiltration with SEC14L2 and SLC6A1 expression in HCC

To explore the role of SEC14L2 and SLC6A1 in immune infiltration in HCC, we evaluated the relationship between differential expression gene and immune cell infiltration by the TIMER platform. A positive correlation between SEC14L2 expression and immune cell infiltration (Cor = 0.172, p = 1.34e-03). The most positive correlation of SEC14L2 and SLC6A1 expression was with B cells (Cor = −0.271, p = 3.33e-07), CD8+ T cells (Cor = −0.224, p = 3.00e-05), macrophages (Cor = −0.325, p = 7.55e-10), neutrophils (Cor = −0.224, p = 2.64e-05), and dendritic cells (Cor = −0.237, p = 1.02e-05) (Figure 8G). We researched the correlation between SEC14L2 and SLC6A1 expression and gene biomarkers of immune cells. The results were most strongly correlated with T cells and T cell exhaustion (Table S1). The correlation between SEC14L2 and gene biomarkers of immune cells in tumor and normal samples is shown in Table S2.

Correlation between methylation and SEC14L2 expression in HCC

The UALCAN database showed that SEC14L2 is highly methylated in the HCC tissue database (Figure 9A). We analyzed the relationship between SEC14L2 methylation and clinical information using the MEXPRESS database. We found significant methylation of SEC14L2 in various clinical factors, including new tumor events after initial treatment, histological type, gender, tumor stage, and OS. The methylation of SEC14L2 occurred on multiple sites, including cg03673688, cg22352499, and cg23665603 (r = 0.441, 0.447, and 0.389, respectively) (Figure 9B). We described the relationship between the methylation sites (cg03673688, cg22352499, and cg23665603) of SEC14L2 and clinical factors of patients using the MethSurv database (Figure 9C).

Figure 9

Methylation of SEC14L2

(A) Methylation of SEC14L2 in HCC tissues was described using the UALCAN database. (B) Relationship between methylation of SEC14L2 and clinical factors by the MEXPRESS database. (C) Methylation sites (cg03673688, cg22352499, and cg23665603) of SEC14L2 described by the MethSurv database.

Construction of the lncRNA-miRNA-mRNA network

To verify the ceRNA mechanism in HCC, we constructed a lncRNA-miRNA-mRNA network after analysis, which included four lncRNAs, one miRNA, and two mRNAs (Figure 10A). We analyzed the correlation between Myc and these key RNAs (Figure S8A). We found that miR-212-3p is a connecting linker gene that participates in ceRNA pathways, whereas lncRNAs indirectly regulated mRNA expression by preferentially encompassing the miRNA response region. We analyzed LINC02691, LINC02499, LINC01354, and NAV2-AS4 and found that lncRNAs were positively correlated with SEC14L2 and SLC6A1 through the same miRNA, miR-212-3p (Figure 10B). Spearman correlation analysis showed that miR-212-3p was negatively regulated in lncRNAs and mRNAs of HCC (Figure 10C). The correlation among lncRNAs, miRNA, and mRNAs is shown (Figure S8B).

Figure 10

Correlation of linear regression analysis between DE lncRNAs and mRNAs

(A) The identified lncRNA-miRNA-mRNA axis is integrated into a circuit map. (B) Relationship of lncRNAs and mRNAs. (C) Correlation between hsa-miR-212-3p and lncRNAs and mRNAs.

Relationship between key RNAs and clinical features

We analyzed the relationship between key RNAs (lncRNAs, miRNAs, and mRNAs) and clinical information, including age, gender, body mass index (BMI), race, TNM stage, distant metastasis, lymph node metastasis, diameter, and prior malignancy (Table S3). SEC14L2 and SLC6A1 had significant differences regarding BMI, tumor stage, and diameter (p < 0.05). TNM, diameter, lymph node metastasis, distant metastasis, and prior malignancy were significantly related to HCC prognosis (p < 0.05). LINC02499 was the lncRNA most significantly correlated with the clinical factors (Table S4).

Discussion

HCC is the fifth-most common cancer, with the second-highest mortality rate worldwide.50, 51, 52 The incidences of liver cancer and mortality have been increasing in the past few years, particularly in Asia. Traditional surgery is no longer the preferred treatment for HCC, so new therapeutic directions need to be explored.,53, 54, 55 Molecular targeted therapy is an important topic for the treatment of HCC; however, the precise molecular mechanisms remain unknown. The development of high-throughput transcriptome sequencing technologies and the ceRNA hypothesis have been proposed, in which lncRNAs and mRNAs interact with each other through shared miRNAs.56, 57, 58, 59 Recently, studies have been conducted to construct a ceRNA network between HCC tissues and adjacent nontumor liver tissues. The molecular mechanisms of ceRNA associated with Myc remain unclear in HCC.

Myc is an oncogene that is dysregulated in >50% of tumors.61, 62, 63, 64 In this study, we first analyzed the large cohort of associated Myc transcriptome profiling with HCC using TCGA database. We selected DE lncRNAs, miRNAs, and mRNAs by comparing the Mychigh tumor tissues with Myclow tumor tissues based on Myc expression, and we constructed a ceRNA network that identified 1125 DE mRNAs, 589 DE lncRNAs, and 93 DE miRNAs.

Numerous studies have reported that lncRNAs play a critical role in different cancers.65, 66, 67, 68 lncRNA MIR503HG inhibits cell proliferation and promotes apoptosis in triple-negative breast cancer (TNBC) cells via the miR-224-5p/HOXA9 axis. lncITPF (long non-coding idiopathic pulmonary fibrosis) promotes pulmonary fibrosis by targeting hnRNP-L depending on its host gene ITGBL1. lncRNA LINC00858 functions through the miR-153-3p/Rabl3 axis, which promotes cell proliferation and infiltration in HCC. Therefore, we analyzed 19 DE lncRNAs in the ceRNA network by multidirectional analysis and calculated risk scores through univariate and multivariate Cox regression analyses. Our results showed that LINC02691, LINC02499, LINC01354, and NAV2-AS4 were significant in predicting OS in HCC, especially LINC02691 and LINC02499. Therefore, the prognostic signature of lncRNAs might act as a prognostic and diagnostic marker for HCC.

Among the prognostic signature lncRNAs, the oncogenic function of LINC01354 has been studied in a variety of cancers—in particular, in colorectal carcinoma, gastric carcinoma, lung carcinoma, and neck squamous cell carcinoma—in which LINC01354 accelerates the proliferation response, migration, and infiltration of cancer cells.72, 73, 74, 75, 76 This study found suppression of LINC01354 expression in Mychigh tumor tissues compared with Myclow tumor tissues. LINC02499 is highly expressed in HCC and inhibits cell-proliferation capacity, migration, and immune cell infiltration in HCC and was first analyzed in the associated Myc ceRNA network. Therefore, we successfully analyzed a series of lncRNAs with functions in the associated Myc ceRNA network. lncRNAs, including LINC02691 and LINC01354, were found to be strongly negatively regulated in HCC through the miR-212-3p/SEC14L2 axis.

miRNA connects lncRNA and mRNA in various cancers, including HCC, acting as a bridge for the ceRNA network.78, 79, 80, 81 miR-96 reduces HCC migration, functioning as a therapeutic target in this disease. miR-221-3p promotes HCC by downregulating the expression of O6-methylguanine-DNA methyltransferase. In this study, we analyzed five miRNAs in the ceRNA network, and only miR-212-3p overexpression was correlated with poor OS. The ceRNA network illustrates that hsa-miR-212-3p may promote cancer cell migration and invasion.84, 85, 86, 87

In this comprehensive study, we screened a ceRNA axis, including the downstream coding gene of SEC14L2 and SLC6A1. We analyzed downstream genes through GO enrichment analysis and KEGG pathway analysis. We found that the GO terms of the dysregulated mRNAs (SEC14L2 and SLC6A1) in HCC could be classified into MF, CC, and BP to further explain the pathways involved. In terms of MF, organic anion transmembrane transporter activity and vitamin E binding suggested that HCC may be a multigene-related disease. During this study, the PPAR signaling pathway was the most important KEGG pathway in the Metascape database, which also has been found in several other cancer types.88, 89, 90, 91, 92, 93, 94

SEC14L2-like phosphatidylinositol transfer proteins SEC14L3/SEC14L2 mediated Wnt/Ca2+ signaling by acting as GTPase proteins. SLC6A1 is overexpressed in prostate cancer and is associated with drug resistance and a poor prognosis but is downregulated in HCC. SEC14L2 and SLC6A1 mRNAs serve as novel transcriptional targets related to Myc and suppress HCC progression and thus may be a promising target for the treatment of Myc-driven HCC. Our results show that Kaplan-Meier survival analysis of ceRNA-correlated genes demonstrated that 19 of the 72 mRNAs had a statistically significant impact on prognosis. Many genes significantly influenced the OS of HCC patients; however, only SEC14L2 and SLC6A1 participated in establishing a complete endogenous regulation network.

This study had some limitations. Because of the lack of other HCC-associated samples, we did not perform further in vitro and in vivo experiments on clinical samples. Additionally, some exploratory experiments remain necessary to identify the functions of the unreported RNAs (lncRNAs, miRNAs, and mRNAs) in this study.

In summary, we introduced a Myc-ceRNA network from genome-wide transcriptome data by various bioinformatics analyses to provide a comprehensive analysis. This approach identified some key RNAs that were significantly associated with prognosis and that provide potential prognostic and diagnostic biomarkers for HCC.

Materials and methods

Data processing and analysis of DE proteins in HCC

We downloaded the transcriptome sequencing data of 421 lncRNAs and mRNAs and 417 miRNAs from 50 patients with HCC through TCGA data portal (https://portal.gdc.cancer.gov/). The 421 RNAs (417 miRNAs) were from 371 tumor tissues and 50 normal tissues. For external validation, we used 10 liver cancer samples and paracancer samples as validation sets from GEO profiles. To study the mechanism underlying the pathogenesis of liver cancer, we defined tumor samples as Mychigh tumor and Myclow tumor with standard median expression levels of MYC. We obtained DE mRNAs, lncRNAs, and miRNAs by comparing the Mychigh and Myclow tumor groups, with thresholds of p < 0.05 and a FC of ≥1.5.

Functional enrichment analysis of DE mRNA

Metascape (http://metascape.org) is an informatics functional annotation tool that integrates many dominating databases to comprehensively analyze genes. We analyzed the enriched function of DE mRNAs using the GO and KEGG pathways in the Metascape website, with the restrictions of p < 0.01, a minimum count of 3, and an enrichment factor of >1.5. We visualized enriched GO pathways by Metascape and presented KEGG analyses using the “ggplot2” package of the R platform.

Cytoscape is a visualized software for network data, which provides users with a more picturesque biological process network. We used lncRNA as a true interaction miRNA target explored by miRcode (http://www.mircode.org/) and used miRDB (http://mirdb.org/) and TargetScan (http://www.targetscan.org/) for miRNA-mRNA target gene prediction. We employed these predicted relationship pairs to construct an endogenous competitive network based on the predicted miRNA expression relationship, which was visualized by Cytoscape version 3.7.0 software. We calculated the densely connected degree of gene nodes in a ceRNA network, and we employed the top 30 genes with the highest confidence scores as hub genes.

Prediction of the subcellular location of lncRNAs

The subcellular localization of lncRNA has an impact on its function. Therefore, we predicted the sequences of the significant lncRNA biomarkers obtained from the LNCipedia (https://lncipedia.org/). Then we used the sequences of lncRNAs to predict the subcellular localization in lncLocator (http://www.csbio.sjtu.edu.cn/bioinf/lncLocator/). A score for each potential subcellular localization of lncRNA included the cytoplasm, nucleus, ribosome, cytosol, and exosome. We analyzed the final results using GraphPad Prism 8.3.

cBioPortal database analysis

cBioPortal (http://www.cbioportal.org/) is an online source application, which provides information on somatic mutations, altered copy number, and mRNA expression, used to visualize cancer genomics data.101, 102, 103, 104 In this study, we used cBioPortal to show SEC14L2 and SLC6A1 genetic changes in HCC.

Immune infiltrate analysis of SEC14L2 and SLC6A1

The TIMER database (http://cistrome.org/TIMER/) is a compositive resource that provides systematic analysis of the abundance of immune cells in the infiltrate of various cancers. We analyzed the expression of SEC14L2 and SLC6A1 in a variety of cancers, as well the correlation between their expression levels and the abundance of the immune cells. In addition, we evaluated the correlation of SEC14L2 and SLC6A1 expression with the biomarkers of immune cells. Gene Expression Profiling Interactive Analysis (http://gepia.cancer-pku.cn/) is an open database that contains tumor and normal samples from TCGA and Genotype-Tissue Expression databases. We used this database to confirm the correlation between DE mRNA and markers of immune cells in tumor and normal samples.

Methylation analysis of SEC14L2

The UALCAN database (http://ualcan.path.uab.edu/) is a data-mining platform, in which the methylation of DE mRNA in tumors can be queried.106, 107, 108 In this study, we used the UALCAN database to analyze SEC14L2 methylation in liver cancer tissues and paracancerous tissues. The MethSurv database (https://biit.cs.ut.ee/methsurv/) is an open website used to obtain CpG methylation data, and it contains significant information about a single CpG. We screened DE mRNA using the MethSurv database and then verified the most important methylated site associated with HCC patient outcomes. MEXPRESS (https://mexpress.be/) is a data-visualization tool used to visualize TCGA expression and the relationship between methylation expression and clinical information.110, 111, 112

The HPA database

The HPA database (https://www.proteinatlas.org/) was initiated in 2003 to map all human proteins in cells, tissues, and organs. In this study, we assembled the protein and RNA expression levels of SEC14L2 and SLC6A1 in various cancers tissues using this database.113, 114, 115

Cancer Cell Line Encyclopedia (CCLE) database analysis

CCLE (https://portals.broadinstitute.org/ccle) is an open-source website, which consists of a large amount of human cancer cell lines, including genomic data, gene expression, copy number, and abundant sequencing data. We analyzed the expression of SEC14L2 and SLC6A1 in HCC cell lines using the CCLE database.

We selected key RNAs according to the previously described analyses, which are presented as the median ± standard deviation. We conducted nonparametric tests to determine whether the expression of RNAs was correlated with the following clinical features: age (≥60 years versus <60 years), gender (female versus male), tumor stage (I + II versus III + IV); BMI (≤18.5, 18.6–23.9, 24–27.9, and ≥28), race (Asian versus White), diameter (≥5 cm versus <5 cm), lymph node and distant metastases (positive versus negative), and prior malignancy (yes versus no). p < 0.05 was used as the cutoff value.

Statistical analysis

We compared differential expression between the Mychigh and Myclow tumor groups using the nonparametric test. The correlation and survival analysis of the relative expression of the important lncRNA-miRNA-mRNA network were processed using GraphPad Prism 8.3. We explored potential lncRNAs using Cox regression analyses in the R package version 4.0. We completed all differential expression analyses using nonparametric tests. p < 0.05 was considered statistically significant.