First of all, the gene transcriptome data (n = 424), gene mutation data (n = 364) and clinical data (n = 377) of patients with HCC were downloaded from the TCGA database (https://portal.gdc.cancer.gov/). Fragments per Kilobase million was used for the transcriptome data, and the VarScan version was used for the gene mutation data. Then, we sorted out the gene expression levels of all lncRNAs from the transcriptome data (normal: 50; tumor: 374) and eliminated the normal samples. Next, we collected survival information for HCC patients from clinical data, excluding samples with incomplete survival information (n = 1) and survival time < 30 days (n = 27) to reduce the interference of non-cancer deaths. After matching the data that met the above conditions, we finally obtained 343 cases with the corresponding gene transcriptome data and clinical data.

Consistent with previous studies (26), we assigned the top 25% and bottom 25% samples in descending order of number of mutations to the Genomically Unstable (GU) and Genomically Stable (GS) groups, separately. Then, the IDs of the samples of the two groups were matched with the LncRNA expression data to obtain the LncRNA expression levels of the GU and GS groups. Finally, Wilcoxon rank-sum test was used on the gene expression levels of the GU and GS groups, and the lncRNAs with a log fold change (logFC) > 1 and a false discovery rate (FDR) < 0.05 were defined as GILncRNAs.

To verify the correlation between the GILncRNAs and GI, we extracted the LncRNA expression data of all tumor samples (n = 374) and conducted a hierarchical cluster analysis of them by calculating the Euclidean distances and cutting the tree into two clusters. We defined one which had a high number of mutations as genomic unstable (GU-like) cluster and the other that had a low number of mutations as genomic stable (GS-like) cluster. Next, we verified the differences between the GU-like and GS-like clusters regarding the number of somatic mutations and the NOD2 ( 11) gene expression levels using the Mann-Whitney U test.

The GILncRNA samples matched with the survival data and expression data were randomly and equally allocated to the training and testing sets. We performed univariate Cox regression analysis on GILncRNAs in the training set and identified GILncRNAs associated with prognosis. Then these GIlncRNAs were analyzed by multiple Cox regression to determine their respective regression coefficients, and a GILncSig was constructed. The formula of the GILncSig was the following:

R i s k s c o r e = ∑ i = 1 n exp ( X i ) × c o e f ( X i ) , where “exp(X_i)”, “coef(X_i),” and “n” represented the expression level, the coefficient, and the GILncSig number, respectively.

Furthermore, we calculated the risk scores of all the samples according to the above formula. The samples with higher or lower average risk scores were assigned to high- and low-risk groups, respectively.

First, we compared the overall survival (OS) and median OS of patients with HCC in the high- and low-risk groups using the log-rank test. We also drew the ROC curves of the 1-, 3-, and 5-year survival rates of the patients with HCC of the testing set and the TCGA set through time-dependent ROC analysis.

Then, we verified the association of the two risk groups with the number of somatic mutations and NOD2 expression levels using the Mann-Whitney U test. Since the genes had wild and mutant states, we analyzed the proportion of the TP53 mutation states between the high- and the low-risk groups using a Chi-squared test. We also implemented a joint survival analysis based on the TP53 mutation status and two clusters using the log-rank test.

Next, we completed univariate and multivariate Cox regression analyses of the risk score and other clinical factors to verify the independence of GILncSig.

Finally, we compared the area under the ROC curve of GILncSig with others for predicting OS at 1-, 2-, and 3-year.

RNAseq data (grade 3) and corresponding clinical information were obtained from the ICGC database (https://dcc.icgc.org/releases/current/Projects) for 240 cases of primary liver cancer (LIRI-JP). Log-rank was used to test the survival difference between the high and low gene expression groups in the KM survival analysis.

A total of 37 cancer and paracancerous tissues diagnosed with HCC were collected from the Second Affiliated Hospital of Nanchang University. This study passed the review of the Medical Ethics Committee of the Second Affiliated Hospital of Nanchang University. The patients included in this experiment were informed and obtained written consent, and this study complied with the standards set by the Declaration of Helsinki.

The human HCC cell line (MHCC97H) was purchased from the Shanghai Institute of Cell Research, Chinese Academy of Sciences, and the cells were validated by the cell bank of short tandem repeats. MHCC97H was cultured in DMEM medium with 10% fetal bovine serum and cultured at 37°C with 5% carbon dioxide. LUCAT1 interfering fragment siRNA and negative control si-NC were obtained from GenePharma (GenePharma Co., Ltd., Shanghai, China). Cell transfection experiments were performed according to the instructions for using Lipofectamine 3000 Transfection Reagent (Invitrogen, Waltham, Massachusetts, USA). The sequences of siRNA are shown in Table S1 .

Total RNA was extracted from tissues and cells using the Trizol method, respectively, and then reverse transcribed into cDNA (TaKaRa, RR047A) and used for real-time quantitative PCR (TAKARA, RR420A). Data analysis was performed using the 2^-ΔΔCt method. The primer sequences used are listed in Table S1 .

The proliferative capacity of MHCC97H with/without LUCAT1 downregulation was observed by CCK-8 assay and EdU assay. Two group cells were seeded in 96-well plates at 6 × 10³ cells per well and allowed to adhere. CCK-8 cell proliferation reagent (10 μl) was added to each well at 0h, 24h, 48h, and 72h. 2 hours after CCK-8 administration, the cell proliferation ability was detected with absorbance detection at 450 nm. In addition, two groups of cells were cultured in a 50 μM EdU medium (Guangzhou Ribobio Co., Ltd.) for 2 h. The medium was then discarded and cells were fixed with 4% paraformaldehyde and 0.5% TritonX-100 in PBS to increase cell membrane permeability. Incubate cells with 1 × Apollo Stain Reaction Solution for 30 minutes at room temperature in the dark. Finally, nuclear DNA was stained with 1 × Hoechst33342 for 30 minutes in the dark. Subsequently, the cells were observed under a fluorescence microscope and photographed.

MHCC97H cells with/without LUCAT1 downregulation were seeded in transwell chambers at 4 × 104 cells per group. Serum-free medium was used for the upper layer of the chamber, and a complete medium with 15% FBS was used for the lower layer of the chamber. After 24 hours of culture, cells in the lower chamber were fixed with 4% formaldehyde solution and stained with 0.2% crystal violet solution for 20 min. Finally, the cells in the upper layer of the chamber were wiped off with a cotton swab, and the cells migrated to the lower layer of the chamber were observed and counted by an inverted microscope. At the same time, cell migration ability was examined using a wound-healing assay. Briefly, two groups of cells were seeded in 6-well plates and scratched with a sterile pipette (200 μL) when the cells were about 90% confluent. The wound healing was then observed with an inverted microscope at 0 hours, 24 hours, and 48 hours, and photographs were taken.

As previously described (27). Total protein from MHCC97H cells transfected with or without si-LUCAT1 was extracted, and western blotting was performed using the following primary antibodies: anti-PDL1/CD274(1:2000, Proteintech, Wuhan, China) and anti-GAPDH(1:10000, Proteintech, Wuhan, China).

The top 10 target genes (mRNAs) that were the most related to GILncRNA were obtained by computing the Pearson correlation coefficients, and the lncRNA-mRNA co-expression network was plotted using the Cytoscape software. We used GSEA (28) to look for the functions and signaling pathways associated with GILncSig. For each analysis, gene set permutations were implemented 1,000 times. A p-value < 0.05 was regarded as the cut-off criteria for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment sorting in GSEA.

The proportions of components in the tumor microenvironment (TME) were assessed using the ESTIMATE algorithm of the “estimate” package (29), resulting in three scores: immune score, stromal score, and ESTIMATE score (30). The higher the score, the greater the proportion of the corresponding component in the TME. We identified potential immune checkpoint genes based on previously published literature (31, 32). The infiltration levels of immune cell types were quantified by ssGSEA in R package gsva (33, 34). ssGSEA is a popular enrichment algorithm, which was extensively utilized in medical studies (35, 36).

To predict chemotherapeutic sensitivity, R-pack “pRRophetic” was used to measure the 50% maximum inhibitory concentration (IC50) of different groups of samples by ridge regression. The IC50 of different groups was then compared by Wilcoxon sign-rank test. Second, immunotherapy responsiveness was predicted using the Tumor Immune Dysfunction and Exclusion (TIDE) Tool (http://tide.dfci.harvard.edu/). As previously shown (37), the GSE78220 dataset was used to analyze efficacy among signature genes and response to immune checkpoint inhibitor treatment.

Wilcoxon rank-sum test was used for comparing the gene expression levels between the GU and GS groups, as well as different tissue samples. Mann-Whitney U test was used to compare the number of somatic mutations and the NOD2 gene expression levels between the two groups. A Chi-square test was used to analyze the clinicopathological characteristics between the training and the testing sets. Log-rank test was used for comparing the OS and median OS. Pearson’s correlation test was used to assess the correlation between GILncSig and immune checkpoint expression levels. p < 0.05 represented a statistical difference.

The general process of this study is shown in Figure 1 . We collated the somatic mutation frequencies for each patient and arranged them in descending order. The patients with the highest and lowest 25% mutation frequencies were then assigned to GU (n = 93) and GS (n = 90) groups, respectively (26). Finally, 88 out of 1,362 lncRNAs were screened and identified as genomic instability-related lncRNAs (GILncRNAs) through the analysis of the gene expression differences between the GU and GS groups ( Table S2 ), which was consistent with previous studies (38). In the GU group, 56 genes showed an increased expression and 32 showed a decreased expression. The 10 most upregulated GILncRNAs and the 10 most downregulated GILncRNAs were observed in the heat maps ( Figure 2A ).

We conducted relevant validation to confirm the correlation between these 88 lncRNAs and GI. In the beginning, we classified all patients into the genomic stable-like (GS-like, n = 163) and genomic unstable-like (GU-like, n = 211) clusters ( Figure 2B ) by unsupervised cluster typing. Then, we found significantly higher somatic mutation frequency and NOD2 (11) gene expression in GU-like clusters compared with GS-like clusters (all p < 0.05, Figures 2C, D ). These results confirmed that the lncRNAs we filtered were correlated with GI.

A total of 343 patients were randomly divided into a training set (n = 172) and a test set (n = 171) in a 1:1 ratio. We screened 13 prognostic-related GILncRNAs from the training set using univariate Cox regression analysis ( Figure 3A ). Then, we performed multivariate Cox regression analysis on the above 13 GILncRNAs, obtained the correlation coefficient of five GILncRNA, and finally established a GILncSig ( Table 1 ). This GILncSig consists of five lncRNAs, including AC245041.2, AP003119.1, MIR210HG, LINC00221, and LUCAT1. The formula of the signature is the following: risk score = expAC245041.2 × 0.086 + expAP003119.1 × 0.173 + expMIR210HG × 0.098 + expLINC00221 × 0.112 + expLUCAT1 × 0.219.

Before checking the signature, we found that the two groups of patients were comparable by comparing the clinical information in the training and testing sets (all p >0.05, Table 2 ).

Next, we tested GILncSig for predictive performance. By plotting the survival curves and ROC curves we found that in the training set, the overall survival(OS) of patients in the high-risk group was significantly lower than that of patients in the low-risk group (median OS = 3.8 versus 8.6 years, respectively, p < 0.001), and the AUCs for predicting patient OS at 1, 3 and 5 years were 0.717, 0.667 and 0.689 ( Figure 3B ) respectively. This was also validated in the testing and TCGA sets ( Figures 3C, D ). In addition, we also plotted the survival and ROC curves of disease-free survival (DFS), disease-specific survival (DSS), and progression-free survival (PFS) of patients in the TCGA set, respectively. The results all showed that the high-risk group had a worse prognosis than the low-risk group, and GILncSig has a good prediction performance (all p <0.05, all AUCs >0.5, Figures 3E–G ).

Furthermore, mutation correlation analysis of the signature was conducted. By drawing a heat map, we found that the GILncSig expression value was upregulated in the high-risk group of the three sets ( Figures 4A–C ). In the mutation point plots and gene expression plots, the number of somatic mutations and the NOD2 expression level showed an increasing trend in the three sets as the risk value increased ( Figures 4A–C ). Through further validation, we found that the number of mutations was significantly elevated in the high-risk groups in the three sets compared with the low-risk groups (all p < 0.05, Figures 4D, F, H ). Similarly, in these groups, NOD2 gene expression levels were higher in the high-risk groups ( Figures 4E, G, I ). In addition, LncRNA has a wild and mutation state; thus, we used the signature to verify the proportion of the single-gene mutation states. In this study, TP53 (39, 40) with a high mutation frequency was selected as the validation object. In the training set, we found that the proportion of TP53 mutation status was significantly higher in the high-risk group compared with the low-risk group(p = 0.002, Figure 4J ). This result was verified in the testing and TCGA sets (all p < 0.001, Figure 4J ). Considering that the different clusters (Gs-like and GU-like) of patients might have different survival rates under different statuses, further analysis was carried out by plotting survival curves. As shown in Figure 4K , for the GU group, there was little difference in survival between TP53 mutant and wild-type patients. In the GS group, we found that the prognosis of patients with TP53 mutation was significantly worse than that of patients in the TP53 wild-type group(p = 0.006). Taken together, these results suggest that our signature can be a good predictor of the frequency of somatic mutations and the prognosis of patients with HCC.

Then, we performed an independent prognostic analysis of GILncSig. Univariate and multivariate Cox regression analysis showed that risk score and tumor stage were independent predictors of OS, DSS, and DFS ( Figures 5A–F ).

At last, to understand the applicability of the GILncSig, we performed a clinical stratification analysis of clinical factors including age and tumor grade. We observed that for patients aged <= 65 or > 65 years with HCC, the survival rate was significantly higher in the low-risk group than in the high-risk group ( Figure 6A , all p < 0.05). Similarly, for patients with tumor grades 1-2 (p = 0.003) or 3-4 (p = 0.051), the survival was higher in the low-risk group than in the high-risk group ( Figure 6B ). This suggests that GILncSig may be generally applicable to patients with HCC.

To more intuitively validate the GILncSig, we compared it with some existing prognostic models for patients with HCC. The first one was the 4-methylated differentially expressed lncRNAs (4-MDELs) model (AC025016.1, LINC01164, LINC01183, LINC01269) established by Liao (41). The second was the 12-lncRNA model constructed by Li (42) (LINC01018, FAM230B, RP11‐290F5.1, APCDD1L‐AS1, FAM99A, LINC01370, AC002116.8, GRM8, RP11‐30J20.1, RP11‐136K7.2, WDR86AS1, and LINC00671). The last was an immune-related LncRNA model created by Cai (43) (FAM120AOS, AL445524.1, AC073257.2, LINC00513, and STK3). The results showed that the AUC values of the ROC of our signature for the 1-, 2-, and 3-year OS were 0.710, 0.653, and 0.634, respectively, which were significantly higher than those of the other three models ( Figures 6C–E ). These results confirmed that our signature was superior to some existing models in predicting the survival rates of patients with HCC. In addition, we confirmed in the ICGC cohort that high expression of the MIR210HG gene in GILncSig could predict poor prognosis in HCC patients ( Figure 6F , p = 0.047).

In GILncSig, LUCAT1 was the most important lncRNA for predicting poor prognosis based on regression coefficients. Therefore, we further evaluated the function of LUCAT1 in HCC. First, expression differential and survival analyses were performed using the ENCORI Online Tools (https://starbase.sysu.edu.cn). It was observed that LUCAT1 was significantly highly expressed in HCC ( Figure 7A , p < 0.001). Patients with HCC in the low LUCAT1 group had a better prognosis compared to the high LUCAT1 group ( Figure 7B , p = 0.012). The qRT-PCR results showed that the mRNA of LUCAT1 gene was highly expressed in HCC tissues compared to the corresponding paraneoplastic tissues ( Figure 7C , n=37, p = 0.002). Further, to analyze the role of LUCAT1, we transfected si-LUCAT1 into MHCC97H cells, the knockdown effect was shown in Figure 7D , and then we used the best knockdown effect si-LUCAT1#2 for the next experiment. To assess the effect of LUCAT1 on proliferation in HCC, we used CCK-8 and EdU staining assays in MHCC97H with/without LUCAT1 knockdown. After interfering with LUCAT1 expression in MHC97H cells, the cell proliferation rate in the si-LUCAT1 group was significantly lower than that in the si-NC group ( Figures 7E, F ). The effect of inhibiting LUCAT1 expression on HCC cell migration was further analyzed. Transwell assays and wound-healing experiments together showed that the migratory ability of MHCC97H cells was significantly reduced upon inhibition of LUCAT1 expression ( Figures 7G, H ). In summary, in vitro experimental data suggest that high expression of LUCAT1 is closely related to poor prognosis in patients with HCC.

To understand the underlying mechanisms via which our signature affected GI, we first mapped the lncRNA-mRNA co-expression network ( Figure 8A ). Then, the possible functions and pathways were discovered using GO and KEGG enrichment analyses. The cellular component (CC) of GO enrichment analysis revealed that the GILncSig-related mRNA was closely connected with the occurrence and development of GI, including the respiratory chain, respiratory chain complex, microtubule-organizing center, attachment site, and meiotic nuclear membrane microtubule tethering complex, among others ( Figure 8B ). Through KEGG pathway analysis, we found that chemical carcinogens-reactive oxygen species, oxidative phosphorylation, HIF-1 signaling pathway, and chemical carcinogenesis, DNA adducts, pentose phosphate pathway, Fructose and mannose metabolism, and other significantly enriched pathways were also significantly associated with GI related. ( Figure 8C ). These results indicated that our GILncSig was involved in GI and its altered expression might destruct the GS of cells and increase GI by interfering with DNA meiosis, microtubule tissue center attachment, and DNA adsorption to affect the normal gene damage repair pathways.

Next, we analyzed the association of GILncSig with the immune status of patients in the TCGA-HCC cohort. High tumor purity has been reported to be associated with immune cell infiltration and cancer development and prognosis (44). Hence, we predicted the stromal score, immune score, and ESTIMATE score in the high- and low-risk groups by the ESTIMATE algorithm. The results demonstrated no difference in immune scores between groups ( Figure 9C , p = 0.76), while there were significant differences in stromal and ESTIMATE scores ( Figures 9A, B , all p < 0.05). Then, we analyzed the correlation between GILncSig and 23 immune cells using the ssGSEA algorithm. The results showed that activated B cells, activated CD8 T cells, CD56 natural killer cells, eosinophils, mast cells, natural killer cells, neutrophils, and type 1 T helper cells were more predominant in the low-risk group, while activated CD4 T cells were more predominant in the high-risk group ( Figure 9D ). It has also been reported that an extremely promising approach to achieve anti-cancer immunity is to block the immune checkpoint pathway (45). Therefore, we analyzed the differences in checkpoint gene expression between the high- and low-risk groups. The result was shown in Figure 9E , which indicated that the expression of all checkpoint genes was statistically significant between the two groups, especially for CD274(PD-L1), CD80, TNFRSF14, HHLA2, HAVCR2, CTLA4, TNFRDF9, TNFRSF9, and CD40 (p < 0.001). These results suggest that GILncSig can be used to assess the tumor immune microenvironment and the expression of immune checkpoint genes in patients with HCC.

To assess the predictive effect of GILncSig on drug therapy for liver cancer, we analyzed the relationship between high- and low-risk groups and the efficacy of commonly used chemotherapeutic agents. Our study showed that the low-risk group was significantly associated with a higher IC50 for chemotherapeutic agents such as bortezomib, gemcitabine, mitomycin, and paclitaxel ( Figures 10A–D , all p<0.001). In contrast, the high-risk group was more sensitive to treatment with axitinib, docetaxel, and lapatinib ( Figures 10E–G , all p<0.01). In immunotherapy, we assessed the effectiveness of immune checkpoint (PD-1 and CTLA-4) inhibitors by TIDE scores in the high- and low-risk groups. TIDE scores were significantly lower in high-risk patients compared to low-risk patients ( Figure 10H , p < 0.001), suggesting that patients in the high-risk group may have better efficacy when receiving immunosuppressive drugs. With the immunotherapy GSE78220 cohort, we also validated that patients with high LINC00221 gene expression in GILncSig responded better to immunosuppressive therapy ( Figure 10I , p=0.0076). Interestingly, we found that inhibiting the expression of LUCAT1 reduced the protein expression of PD-L1 by western blotting ( Figure 10J , p = 0.0024).