Comprehensive analysis of an immune-related ceRNA network in identifying a novel lncRNA signature as a prognostic biomarker for hepatocellular carcinoma

DElncRNAs, DEmRNAs, and DEmRNAs in HCC

By comparing HCC and adjacent normal liver tissues in TCGA-LIHC database (p < 0.05, |log₂-fold change (FC)| > 1), a total of 3724 DElncRNAs (3249 upregulated and 475 downregulated), 5199 DEmRNAs (4088 upregulated and 1111 downregulated), and 330 DEmiRNAs (287 upregulated and 43 downregulated) were identified (Figure 2A–2C).

Constructing the lncRNA-miRNA-mRNA ceRNA network of HCC

The lncRNA-miRNA pairs containing 332 lncRNAs (286 up-regulated and 46 down-regulated) and 30 DEmiRNAs were identified in the miRcode online database. The target gene predictions of 30 DEmiRNAs were identified by using following online tools: TargetScan, miRDB and miRTarBase. The interaction pairs of miRNA–mRNA involving 30 DEmiRNAs and 1494 DEmRNAs were confirmed. The predicted target gene was matched with DEmRNAs, and a total of 242 target DEmRNAs were obtained. After the removal of the 1 remaining miRNA-lncRNA pair, 331 DElncRNAs, 29 DEmiRNAs, and 242 DEmRNAs were finally obtained and the ceRNA network of HCC was built based on them (Supplementary Figure 1 and Supplementary Table 1).

Constructing and validating the ceRNA-and-immune-related lncRNA prognostic signature

Pearson correlation analysis revealed that 172 lncRNAs out of the 331 ceRNA-related DElncRNAs showed correlation with 1308 immune-related mRNAs (|R| > 0.4 and P < 0.001). The 172 lncRNAs were identified as ceRNA-and-immune-related lncRNAs to construct an immune-related prognostic signature (Supplementary Table 2). In the training cohort, the 172 lncRNAs were subjected to univariate Cox regression analysis and least absolute shrinkage and selection operator (LASSO) Cox regression analyses with log λ ≈ -3.10. Then 11 prognosis-related lncRNAs were filtered out (Figure 3A–3C and Table 1). The risk score of the prognosis signature was calculated as follows:

Risk score = (0.102670442*expression value of AP002478.1) + (0.055778189*expression value of C10orf91) + (0.035841281*expression value of LINC00501) + (0.008492106*expression value of ERVMER61-1) + (0.026966520*expression value of AC114489.1) + (0.131909917*expression value of AL118511.1) + (0.048942893*expression value of MIR137HG) + (0.059437958*expression value of SERHL) + (0.007189903*expression value of LINC00346) + (0.003028755*expression value of LINC00491) + (0.005353952*expression value of DSCR4-IT1).

The optimal cutoff value of the risk score in the training cohort, which divided 181 HCC samples in the training cohort into the high-risk group (n = 69) and the low-risk group (n = 112), was 1.00 (Figure 3D). According to the results which were illustrated in Figure 4A, 4B, higher risk scores were associated with more deaths (Figure 4A, 4B). And the 11 lncRNAs expression levels were upregulated in the high-risk group (Figure 4C). According to the Kaplan-Meier survival curve in Figure 4D, the low-risk group had longer survival time than the high-risk group (P < 0.0001). The receiver operating characteristic (ROC) curve suggested that the risk score could predict the 1-year (0.838), 3-year (0.72) and 5-year (0.736) survival rate in HCC patients effectively (Figure 4E). Furthermore, as illustrated in Figure 4F, multivariate Cox regression analyses suggested that the risk score of the signature (HR = 1.951, P < 0.001) could be the independent prognostic indicator of overall survival (OS).

The survival model with the same cutoff value (1.00) of the training cohort was applied to the testing cohort (n = 187) and the entire TCGA-LIHC cohort (n = 368), and the validation results of the signature in these two cohorts also demonstrated excellent capacity for survival prediction (Supplementary Figures 2, 3).

Correlation of clinical characteristics with the immune-related lncRNA prognostic signature

To assess the wide applicability of the signature, the entire TCGA-LIHC cohort was classified into different groups by different clinical characteristics to perform the correlation analysis. The stratification analysis results suggested that the clinical characteristics exhibited distribution patterns were consistent with the risk score of the signature (Figure 5A). The ROC curves for 5-year OS manifested that the risk score (AUC = 0.725) was better than other clinical characteristics (Figure 5B), indicating the considerable reliability of the immune-related lncRNA signature for HCC.

Cancer hallmark and immune signature analysis

Hallmark enrichment analysis (h.all.v7.2.symbols) in this study revealed that the high-risk group was primarily linked with G2M checkpoint, PI3K/AKT/MTOR signaling pathway, glycolysis, and some canonical signaling pathways (Figure 6A–6C). These findings might help to understand the biological functions of the ceRNA-and-immune-related signature. In addition, our signature also regulated some immunologic signatures (c7.all.v7.0.symbols), such as germinal center (GC) B cell vs. plasma cell up, memory vs. naïve CD8 T-cell down and CD4 T-cell vs. NK T-cell down, demonstrating that the lncRNA signature might play an essential role in immune-related regulation (Figure 6D–6F).

Relation of the lncRNA signature with immune cell infiltration and ICB

The prognostic lncRNA signature was significantly positively correlated with immune infiltration of eosinophils (R = 0.16; p = 0.0022), M0 macrophage (R = 0.33; p < 0.001), follicular helper T cells (R = 0.12; p = 0.026), and regulatory T cells (R = 0.31; p < 0.001). And the signature was also significantly negatively correlated with immune infiltration of resting mast cell (R = -0.18; p < 0.001), monocytes (R = -0.18; p < 0.001), and CD4 memory resting T cells (R = -0.25; p < 0.001) (Figure 7A–7G). The high correlation between the signature and immune cell infiltration strongly manifested that the signature was consistent with the immunologic signature analysis.

The ICB correlation analysis, including PD-1, CTLA4, TIM-3, PD-L1, PD-L2, IDO1, GITR, SOAT1, CDK1, HDAC2 and MMP9, was conducted with the signature. As shown in Figure 8, the lncRNA signature showed significantly positive correlation with all ICBs targeting molecules. This result suggested that the signature might contribute to the evaluation of response to ICB therapy in HCC.

Identification of the hub gene in the lncRNA signature

LINC00491 has already been recognised as hub oncogene in non-small-cell lung cancer (NSCLC) and CA, but it has never been explored and validated in HCC [9, 10]. Therefore, the role of LINC00491 with regard to HCC requires further assessment.

As shown in the boxplot, LINC00491 expression level was significantly upregulated in the HCC tissues, and LINC00491 also significantly differentially expressed in different groups of varying vital status (Figure 9A, 9B). Interestingly, LINC00491 level had the ability to distinguish noncancerous tissues from HCC tissues at AJCC stage I, II, III, and IV (Figure 9C). And LINC00491 could also be applied to distinguish normal tissues from tumor tissues at different histologic grades (Figure 9D). To explore the prognostic ability of LIHC00491, survival analysis found that HCC patients with higher LINC00491 expression level demonstrated a poorer OS (Figure 9E). Subgroup results showed that patients with higher LINC00491 expression level in histologic stage G3/G4, AJCC stage I/II, male and old-age (age ≥ 60) had shorter OS (Figure 9F–9M). Following the multivariate Cox analysis result, LINC00491 expression level could act as an independent prognostic indicator of OS in HCC patients (Figure 9N). Furthermore, LINC00491 was also positively correlated with the above mentioned 11 ICBs (Supplementary Figure 4).

LINC00491 in vitro and in vivo experiments

LINC00491 was overexpressed in cancerous tissues by qRT-PCR in Zhujiang cohorts, which was consistent with the above database results (Figure 10A). To illuminate the proliferative rate after LINC00491 knockdown in HCC cells, Cell Counting Kit-8 (CCK-8) assays in MHCC-97H and HLF cell lines both showed lower absorbance at 96h and 120h, respectively (Figure 10B). After silencing of LINC00491, the Colony formation assay indicated that the growth of HCC cells was significantly suppressed after silencing of LINC00491 (Figure 10C), transwell assay results showed that the invasion and migration rates of HCC cells were significantly suppressed (Figure 10D, 10E), and wound healing in HCC cells were also significantly inhibited (Figure 10F). Cell cycle progression in MHCC97H and HLF cells revealed that LINC00491 depletion significantly repressed cell cycle progression (Figure 10G). In vivo experiments with subcutaneous tumor models were shown in Figure 10H. In accord with in vitro observations, tumor volume and tumor weight were both significantly decreased after knockdown of LINC00491 (Figure 10I, 10J).

Constructing the immune-related ceRNA network in HCC

242 DEmRNAs in the above-mentioned HCC ceRNA network were intersected with the IMMPORT gene list and 33 immune-related DEmRNAs were obtained. To obtain more reliable results, a correlation analysis (R > 0.4, p < 0.001) was carried out between former 172 ceRNA-and-immune-related DElncRNAs and the 33 immune-related DEmRNAs. The results confirmed 161 pairs of lncRNA-mRNA containing 63 lncRNAs and 26 miRNAs. By matching the relationship of the above-mentioned ceRNA network, the immune-related ceRNA network was finally built by 59 lncRNAs, 21 miRNAs, and 26 mRNAs (Figure 11 and Supplementary Table 3).

Differential expression analyses

TCGA-LIHC data portal as the HCC database contained 374 HCC-tissue samples and 50 adjacent normal liver-tissue samples. EdgeR package was used to identify DElncRNAs, DEmRNAs and DEmiRNAs between HCC tissues and normal tissues [22]. |log₂FC| > 1 and p < 0.05 were selected as significant cutoff values. This study does not require the approval of the ethics committee.

Construction of ceRNA network

The miRcode database was used to predict the relationship between DElncRNAs and DEmiRNAs. The miRDB, miRTarBase, and TargetScan databases were applied to retrieve the target DEmRNAs of DEmiRNAs. Cytoscape software was used to visualize the ceRNA network.

Identifying immune-related lncRNAs and constructing prognostic signature

The Immunology Database and Analysis Portal (ImmPort) database contained a total of 1811 immune-related genes. Pearson correlation analysis using the thresholds of Spearman's correlation coefficient with an absolute value (|R|) of > 0.4 and p < 0.001 was performed to identify immune-related lncRNAs in the above ceRNA network of HCC.

Based on a 5 : 5 ratio, the entire TCGA-LIHC cohort (n = 368) with prognostic clinical information was randomly divided into the training (n = 181) and the testing (n = 187) cohorts. Univariate Cox regression analysis and LASSO Cox regression analysis were utilized in the training cohort to identify prognostic immune-related lncRNAs and develop a prognostic signature [23, 24]. The optimal cutoff value of risk scores calculated by survminer package (http://cran.r-project.org/) in R software was used to classify the training cohort into the low- and the high-risk groups. And to further confirm the predictive accuracy, the signature was also applied to the testing and the entire TCGA-LIHC cohorts.

Exploring the bio-function of the lncRNA signature

To find enriched function terms, gene set enrichment analysis (GSEA) with the threshold of p < 0.05 and false discovery rate (FDR) < 0.25 was performed using the Molecular Signatures Database (MSigDB) [25]. CIBERSORT was used to analyze the relation between the immune cell infiltration and the lncRNA signature in HCC samples [26]. To explore the relation with ICB therapy, the correlation analysis was performed between the immune-related signature and 11 immune checkpoint inhibitors, including PD-1, CTLA4, TIM-3, PD-L1, PD-L2, IDO1, GITR, SOAT1, CDK1, HDAC2 and MMP9 [27–35].

The material preparation of in vitro and in vivo experiments

28 pairs of fresh frozen specimens with corresponding normal tissues detected by qRT-PCR were obtained from HCC patients of the Department of Hepatobiliary Surgery of Zhujiang Hospital. The Medical Ethics Committee of Zhujiang Hospital had approved our study. According to the manufacturer’s protocol, Lipofectamine 3000 (Invitrogen, Carlsbad, CA, USA) was purchased for transient transfection of LINC00491-Smart Silencer (si-LINC00491) and si-normal control (si-NC) (Ribobio, Guangzhou, China) in MHCC-97H and HLF cell lines. And to conduct animal experiments, the specific short hairpin RNAs to LINC00491 (sh-LINC00491) was used to transfect HLF cells.

Migration assay

For the transwell assay, cells suspended in serum-free medium were placed in the upper chamber of 24-well Transwell plates (Corning, NY, USA), and the culture medium containing 10% FBS was seeded into the lower chamber. For wound healing assay, the cells were plated into 6-well-plates and cultured in DMEM without FBS. A 200-μl pipette tip was used to scrape the monolayer. At 0 and 48 hours after scratch generation, cell migration was recorded with white light using a phase contrast microscope.

Proliferation assay, colony formation assay and cell cycle analysis

The CCK-8 assay (Vazyme Biotech, Nanjing, China) was used to determine the proliferative capacity of cells. According to the manufacturer’s instructions, cells were plated in 96-well plates for 4 days, and the absorbance was observed at 450 nm. For colony formation assay, cells were plated in 6-well plate and they were cultured in DMEM containing 10% FBS at 37° C for 2 weeks. Before counting, the colonies were fixed in methanol, and then the colonies were stained with 0.1% crystal violet for 30 min. Cell cycle analysis was carried out as previously described [36].

Mouse xenograft assay

To perform cells transfection, the left and right backs of nude mice (4-week-old male BALB/c mice) were subcutaneously injected with sh-LINC00491 or sh-NC cells, respectively. The growth of tumors was recorded every 3 days. Mice were killed on the 3 weeks after inoculation, and the resected tumors were weighed for subsequent analysis. All animal experiments were approved by the Ethics Committee for Laboratory Animals of Zhujiang Hospital, Southern Medical University, Guangzhou, China.