Identification of an EMT-related lncRNA signature and LINC01116 as an immune-related oncogene in hepatocellular carcinoma

Data collection and processing

Gene expression data of 424 samples including 374 tumor and 50 normal samples and clinical information of 374 HCC patients were downloaded from TCGA database. LncRNAs and mRNAs were identified from RNA sequence data based on annotations in the human GENCODE database. The “edgeR” package in R software was used to identify differentially expressed lncRNAs with the threshold P value <0.05 and |log(fold change)|>2. Additionally, the data of 368 HCC patients with complete clinical information were randomly divided into training and validation sets at a 1:1 ratio. The clinical characteristics of the 2 groups were similar.

EMT-related lncRNA identification

EMT-related genes were obtained from the HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION gene set in Molecular Signatures Database v7.4 (https://www.gsea-msigdb.org/gsea/msigdb/). EMT-related lncRNAs were identified by correlation analysis based on a P value <0.001 and absolute Pearson coefficient >0.4.

Construction of the EMT-related lncRNA signature

A flow diagram of the study is shown in Figure 1. In order to establish an EMT-related lncRNA signature for predicting prognosis in HCC, we carried out a univariate Cox regression analysis to screen prognostic lncRNAs (P<0.05), followed by least absolute shrinkage and selection operator (LASSO) penalized Cox regression and multivariate Cox regression analyses to select those that were relevant. The lncRNA expression level and regression coefficient β from the multivariate Cox regression results were used to develop the following equation: risk score=β lncRNA1×lncRNA1 expression+β lncRNA2×lncRNA2 expression+/+β lncRNA n×lncRNA n expression. After assessing the prognostic value of the signature, the “Survminer” package in R was used to determine the optimal cutoff value for dividing patients into high- and low-risk groups.

Figure 1. The flow chart of this study.

Validation of the prognostic significance of EMT-related lncRNAs

The prognostic value of the EMT-related lncRNA signature was verified by Kaplan–Meier and receiver operating characteristic (ROC) curve analyses in the training and validation sets and in the combined dataset (i.e., the whole patient population). The relationship between risk score and survival is presented as a risk score distribution map and lncRNA expression heatmap. Combined with other clinical characteristics of HCC patients, uni- and multivariate Cox regression analyses were carried out to determine whether the risk score was an independent prognostic factor in HCC. The prognostic accuracy of the signature and other clinical features in all groups was verified by ROC curve analysis.

Analysis of the relationship between the prognostic signature, immune cell infiltration, and immune checkpoint targets

The CIBERSORT algorithm was used to analyze the abundance of 22 types of tumor-infiltrating immune cell (TIIC) and Tumor Immune Estimation Resource (TIMER) was used to calculate the abundance of 6 different types of TIIC. The association between the lncRNA signature and the different types of TIIC was analyzed to determine whether the signature can predict immune cell infiltration in HCC.

The success of programmed death (PD)-1/programmed death ligand (PD-L)1 and cytotoxic T-lymphocyte antigen (CTLA)-4 blockade for HCC treatment has intensified research efforts to identify novel immune checkpoint targets. In this study, we examined established targets as well as those with therapeutic potential such as T-cell immunoglobulin and mucin domain-containing (TIM)-3, indoleamine 2,3-dioxygenase (IDO)1, glucocorticoid-induced tumor necrosis factor receptor-related (GITR), histone deacetylase (HDAC)2, B7-histone (H)3, and V-domain Ig suppressor of T cell activation (VISTA) [13–15]. In order to evaluate the role of the lncRNA signature in the treatment response to immune checkpoint blockade (ICB), we performed a correlation analysis between the signature and the expression level of 10 immune checkpoint targets in HCC.

Tissue specimens and cell culture

We collected 50 paired HCC tissue and adjacent nontumor tissue samples from patients who underwent hepatectomy between March 2016 and October 2017 at the Hepatic Surgery Center of Tongji Hospital. OS information was collected through electronic medical records or telephone follow-up. Informed consent forms to donate their tissue samples for biomedical research, which were approved by the ethics committee of Tongji Hospital, were signed by all patients. Huh7 and SK-Hep1 HCC cell lines were purchased from the Liver Cancer Institute of Fudan University (Shanghai, China). The cells were cultured in Dulbecco’s Modified Eagle’s Medium supplemented with 10% fetal bovine serum at 37° C and 5% CO₂.

Immunofluorescence analysis

HCC cells were cultured on coverslips in 24-well plates for 24 h. After 3 washes with phosphate-buffered saline, the cells were fixed with 4% paraformaldehyde for 15 min and permeabilized with 0.1% Triton X-100 for 10 min. After culturing in medium containing 1% bovine serum albumin for 2 h, the cells were incubated overnight at 4° C with rabbit anti–E-cadherin and mouse anti-vimentin antibodies and then labeled with Alexa Fluor 568- or 488-conjugated secondary antibody (from different species) for 2 h. Cell nuclei were stained with 4′, 6-diamidino-2-phenylindole for 10 min. The coverslips were mounted and observed with a confocal microscope.

Fluorescence in situ hybridization (FISH)

Cells were cultured on coverslips in 24-well plates. LINC01116 was detected in HCC cells by FISH using a probe and commercial kit from RiboBio (Guangzhou, China) according to the manufacturer’s instructions. The labeled cells were observed with a confocal microscope.

Xenograft tumor model

Male BALB/c nude mice (4 weeks old, weighing 15 g) were maintained under specific pathogen-free conditions. Huh7 cells transfected with empty vector or LINC01116 overexpression plasmid (3×10⁶) were subcutaneously injected into the back of the mice on the right and left sides, respectively. Tumor growth was observed every 3 days and tumors were harvested after 3 weeks.

Immunohistochemistry

Tumor tissues from mice were fixed in 4% paraformaldehyde solution and then dehydrated in ethanol solutions, embedded in paraffin, and cut into sections that were incubated overnight at 4° C with primary antibodies against Ki67, vimentin, and E-cadherin followed by horseradish peroxidase-conjugated secondary antibody for 1 h at room temperature. Images of representative fields were acquired with a C1 Digital ECLIPSE microscope.

Cell proliferation and migration assays and cell cycle analysis

The cell proliferation (Cell Counting Kit [CCK]-8), 5-ethynyl-2′-deoxyuridine (EdU), transwell migration, and wound healing assays and cell cycle analysis were carried out as previously described [16].

Statistics analysis

Data analysis were conducted using R software (version 3.6.3) and Prism 7.0 software. Data are presented as mean ± SD. Qualitative variables were analyzed by Student's t test, one-way ANOVA and Pearson correlation test were used to analyze qualitative variables. P<0.05 was considered to be statistically significant.

Data availability statement

Publicly available datasets were analyzed in this study. The data are accessible at the following repositories: https://portal.gdc.cancer.gov/repository.

Ethics statement

This research protocol was approved by the ethical committee of Tongji Hospital and a consent form was signed by each participating patient.

Identification of EMT-related lncRNAs

Using a cutoff of |log₂(fold change)|>2 and P<0.05, we identified 1283 differentially expressed lncRNAs (84 downregulated and 1199 upregulated) in 374 HCC samples compared to 50 normal samples. The top 50 lncRNAs were used to construct a heatmap (Figure 2A). We also identified 200 EMT-related genes from the Molecular Signatures Database; by correlation analysis, 335 EMT-related lncRNAs were identified based on a threshold of P<0.001 and absolute Pearson coefficient >0.4.

Figure 2. (A) The top 50 of both up- and down-regulated DElncRNAs with the most significant differences were represented by the heatmap. (B) EMT-related lncRNA signature was established by Univariate Cox regression, LASSO and Multivariate Cox regression analysis. (C) Kaplan–Meier survival analysis suggested a lower OS in high-risk groups. ROC curve showed good accuracy of this signature in predicting OS of 1, 3 and 5 years.

Construction of the EMT-related lncRNA signature in HCC

A total of 368 TCGA HCC patients with complete clinical information were randomly assigned to the training (n=197) and validation (n=171) sets. In order to identify lncRNAs with prognostic value in HCC, the EMT-related lncRNAs were evaluated by univariate Cox analysis; 180 were found to be significantly related to the overall survival (OS) of HCC patients. LASSO and multivariate Cox regression analyses were performed to select and validate the prognostic lncRNAs in the training cohort (Figure 2B). An EMT-related prognostic lncRNA signature for HCC was established based on the regression coefficients and expression levels of these lncRNAs that included AC011352.3, AC009005.1, AC025176.1, AP002478.1, and LINC01116.

The Kaplan–Meier analysis of the training and validation sets and whole patient population indicated that HCC patients with a high risk score had a shorter OS than those with a low risk score. The ROC curve analysis for the prognostic lncRNA signature showed excellent predictive accuracy, with an area under the ROC curve (AUC) of 0.786, 0.717, and 0.736 for predicting 1-, 3-, and 5-year OS, respectively. The prognostic value of the signature was assessed in the validation set and whole patient population; the AUCs for predicting 1-, 3-, and 5-year OS were 0.759, 0.757, and 0.754, respectively, in the validation set and 0.777, 0.741, and 0.738, respectively, in the whole population (Figure 2C).

Validation of the prognostic lncRNA signature

The mortality rate in all groups increased with risk score; meanwhile, the expression levels of the 5 hub lncRNAs in the EMT-related signature also increased (Figure 3A). A multivariate Cox regression analysis that also included other clinical characteristics of HCC patients showed that risk score was an independent prognostic factor in HCC, with a hazard ratio of 1.690 (1.204−2.372) in the training set, 1.503 (1.159−1.948) in the validation set, and 1.446 (1.229−1.702) in the whole patient population (Figure 3B). The risk score had a higher AUC than the other clinical features in each group, confirming its value as an accurate prognostic factor (Figure 3B).

Figure 3. (A) Mortality status and lncRNAs expression in each patient were plotted according to the ordered risk score. (B) Multivariate and ROC curves confirmed the EMT-related lncRNA signature as an independent prognostic factor for HCC.

To investigate the robustness of the prognostic lncRNA signature, we evaluated its predictive accuracy in the Tongji cohort including 50 HCC patients. Consistent with the results from the TCGA cohort, that HCC patients with a high risk score had a shorter OS than those with a low risk score (Figure 4A). In addition, the prognostic lncRNA signature had an AUC of 0.754 at 1 year, 0.723 at 3 years, and 0.703 at 5 years in the Tongji cohort (Figure 4B), indicating an excellent predictive accuracy. Then, we analyzed the relationships between the EMT-related lncRNA signature and clinicopathologic characteristics of HCC patients. The risk score was significantly higher in patients with vascular invasion and those who died (Figure 4C, 4D). We also observed that the risk score tended to increase with tumor grade and stage (Figure 4E, 4F), implying that the signature reflects the high malignancy of HCC.

Figure 4. (A) Kaplan–Meier survival analysis in Tongji cohort. (B) ROC curve showed good accuracy of this signature in Tongji cohort. (C) The association between risk score and living status. (D) The association between risk score and vascular invasion. (E) The association between risk score and tumor grade. (F) The association between risk score and tumor stage. **P < 0.01 and ***P < 0.001.

Functional analysis of the prognostic signature

We performed a principal component analysis to compare low and high-risk patients in the training and validation sets and the whole patient population who were classified based on the EMT-related lncRNA signature. The low- and high-risk groups formed distinct clusters, confirming the prognostic accuracy of the signature (Figure 5A). We also explored the function of lncRNAs in the signature by performing a correlation analysis with mRNAs that were differentially expressed between high- and low-risk groups, whose functions were predicted by Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway analyses. It showed that function annotations of this signature was enriched in terms related to signal release, intrinsic component of synaptic membrane, hormone activity and neuroactive ligand−receptor interaction (Figure 5B, 5C).

Figure 5. (A) PCA showed a scattered distribution of patients in each group. (B, C) GO and KEGG analysis of DEGs between high and low-risk groups of the entire group.

Correlation between the prognostic signature, immune cell infiltration, and immune checkpoint targets

Using the CIBERSORT algorithm, we investigated whether the EMT-related lncRNA signature was related to immune cell infiltration in HCC. The results of the correlation analysis revealed a significant correlation with infiltrating resting memory cluster of differentiation (CD)4+ T cells (r=−0.154; p=3.08e−03), regulatory T cells (r=0.197; p=1.46e−04), monocytes (r=-0.149; p=4.27e−03), macrophages M0 (r=0.265; p=2.53e−07), resting mast cells (r=−0.163; p=1.69e−03), and eosinophils (r=0.171; p=9.73e−04; Figure 6A). Additionally, all lncRNAs in the signature showed a significant correlation with the infiltration of multiple immune cell types.

Figure 6. (A) Correlation analyses of the EMT-related lncRNA signature with immune cell infiltration. (B) Correlation analyses of the EMT-related lncRNA signature with immune checkpoint targets.

To investigate the potential role of the EMT-related lncRNA signature in ICB therapy for HCC, we performed a correlation analysis between risk score and 10 immune checkpoint targets (CTLA4, TIM-3, PD-L1, PD-L2, PD-1, IDO1, GITR, HDAC2, B7-H3, and VISTA) and found close correlations (Figure 6B), suggesting that the signature can predict the outcome of ICB therapy in HCC patients.

Hub lncRNAs in HCC

Using TCGA data, the expression levels of the 5 hub lncRNA s in this signature were compared in the HCC tissues and normal liver tissues. It showed that all the 5 hub lncRNAs expression levels were higher in HCC tissues than these in the normal liver tissues. Besides, the Kaplan–Meier analysis indicated that the increased the 5 hub lncRNAs were significantly correlated with OS of HCC patients (Supplementary Figure 1). The increased levels of the 5 hub lncRNAs in HCC was confirmed in HCC specimens from the Tongji Hospital cohort, and the Kaplan–Meier analysis with the log-rank test indicated a significant correlation with OS (Figure 7A). Additionally, 3 of the lncRNAs (AC025176.1, LINC01116, and AC011352.3) were highly expressed in 4 HCC cell lines compared to a normal liver cell line (Figure 7B). Of these, LINC01116 was selected for further analysis based on the fact that it had the highest relative expression in clinical HCC specimens and has not been previously characterized in HCC.

Figure 7. (A) The expression levels of 5 hub lncRNAs in HCC specimens of Tongji hospital. (B) The expression levels of 5 hub lncRNAs in normal liver cell line and 4 HCC cell lines.

LINC01116 enhances HCC cell proliferation and regulates the cell cycle

To determine the subcellular localization of LINC01116 in HCC cells, we performed cytoplasmic/nuclear fractionation and FISH in SK-Hep1 and Huh7 cells. LINC01116 was detected at a higher level in the cytoplasm than in the nucleus, implying that it mainly acts on cytoplasmic factors (Figure 8A, 8B).

Figure 8. (A) The expression level of LINC01116 in the subcellular fractions of HCC cells was detected by qRT-PCR. (B) FISH assay analysis for the location of LINC01116 in HCC cells. (C) Transfection efficiency was verified by qRT-PCR. (D–F) Cell viability was evaluated with CCK-8, EdU and colony formation assays in HCC cells. (G) Cell cycle was examined by the flow cytometry. (H) Transwell assays were used to detect HCC cells invasion and migration. (I) Migration ability was evaluated by wound healing assay. **P < 0.01 and ***P < 0.001.

To investigate the function of LINC01116 in HCC, we used 2 different small interfering (si)RNAs to silence its expression in SK-Hep1 cells. Conversely, LINC01116 expression was markedly increased in Huh7 cells transfected with LINC01116 overexpression plasmid (Figure 8C).

The effect of LINC01116 loss or gain of function on HCC cell proliferation and viability was evaluated with the CCK-8, EdU, and colony formation assays. The results showed that LINC01116 knockdown suppressed the proliferation of SK-Hep1 cells while LINC01116 overexpression enhanced Huh7 cell proliferation (Figure 8D–8F). We also examined cell cycle distribution under the 2 conditions by flow cytometry. LINC01116 silencing increased the proportion of SK-Hep1 cells arrested in G0/G1 phase whereas LINC01116 overexpression induced cell cycle progression in Huh7 cells (Figure 8G). In contrast, HCC cell apoptosis was unaffected by LINC01116 depletion and overexpression (Supplementary Figure 2).

LINC01116 promotes EMT and metastasis in HCC

To investigate the role of LINC01116 in the migration and invasion of HCC cells, we performed transwell migration and wound healing assays. LINC01116 knockdown suppressed migration and invasion in SK-Hep1 cells, while its overexpression had the opposite effect in Huh7 cells (Figure 8H). In the wound-healing assay, migration distance was decreased in SK-Hep1 cells transfected with siRNAs and was increased in Huh7 cells overexpressing LINC01116 (Figure 8I). We also examined the expression of the EMT markers vimentin, E-cadherin, and N-cadherin by western blotting and immunofluorescence analysis. LINC01116 knockdown in SK-Hep1 cells decreased the levels of vimentin and N-cadherin and increased that of E-cadherin, with LINC01116 overexpression have the opposite effects in Huh7 cells (Figure 9A, 9B).

Figure 9. (A) EMT markers were examined via western blot analysis. (B) Immunofluorescence assay was used to detect EMT markers. (C) Image of subcutaneous tumor tissues. The volume and weight of tumors were measured. (D) Ki67, vimentin and E-cadherin were observed in subcutaneous tumor tissues by IHC. **P < 0.01 and ***P < 0.001.

To verify the above results in vivo we used a mouse xenograft model. Compared to the control group, Huh7 cells overexpressing LINC01116 formed larger tumors (Figure 9C); moreover, these tumors showed upregulation of Ki67 and vimentin and downregulation of E-cadherin by immunohistochemistry (Figure 9D). Collectively, these results indicate that LINC01116 promotes EMT and metastasis in HCC.

LINC01116 is associated with immune cell infiltration and expression of immune checkpoint targets

As described above, LINC01116 expression was found to be significantly correlated with the infiltration of 6 immune cell types based on the CIBERSORT algorithm (Figure 6A). We further examined the relationship between LINC01116 and TIICs using TIMER. LINC01116 expression level was positively and significantly correlated with B cell, CD4+ T cell, CD8+ T cell, macrophage, neutrophil, and dendritic cell infiltration (Figure 10A). To evaluate the influence of LINC01116 on ICB therapy in HCC, we performed a correlation analysis between LINC01116 and immune checkpoint targets (CTLA4, TIM-3, PD-L1, PD-L2, PD-1, IDO1, GITR, HDAC2, B7-H3, and VISTA) and found positive correlations with the levels of all 10 targets (Figure 10B). Thus, LINC01116 plays a critical role in immune regulation in HCC and is a potential target for ICB therapy.

Figure 10. (A) Correlation analyses of LINC01116 with immune cell infiltration based on TIMER database. (B) Correlation analyses of LINC01116 with immune checkpoint targets.