Both expression data (lncRNA, mRNA) and clinical data of HCC patients were downloaded from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) website. We integrated the downloaded gene expression profile data with perl script and converted Ensembl ID to gene ID to obtain expression matrix and clinical file. The expression level of m6a gene was extracted from the transcriptome matrix, and 23 m6a-related lncRNAs were found by co-expression analysis using the limma package. According to the expression level of m6a-related lncRNA in the sample, the network relationship file and the node attribute file, the igraph package was used to visualize the co-expression network. Then, the expression data of m6a-related lncRNA and the survival data extracted from the clinical data were combined, and the survival package was used for univariate Cox analysis. The obtained prognosis-related lncRNAs were intersected with differentially expressed lncRNAs between tumor and non-tumor tissues to obtain prognosis-related differential lncRNAs.

A total of 370 HCC patients were clustered into different clusters by ConsensusClusterPlus R package and finally determined by consensus cumulative distribution function and delta area. Then, survminer package was used to perform survival analysis on different classification groups according to the classification results and survival data, and the survival curve was obtained. The expression of other clinical parameters and prognosis-related lncRNA were extracted, and the heatmap package was used to draw a heat map of the correlation between classification and clinical. Correlation heatmap between CD274 and prognosis-related lncRNAs was performed after differential analysis of CD274 by limma R package. In addition, heatmap between clinical data, including age, gender, grade, stage, tumor, node, metastasis as well as clustering group, and down-/up-regulation of lncRNAs was constructed as well.

According to the gene expression matrix, the cibersort algorithm was used to analyze the immune cell content file and the difference analysis was performed. The filter condition p value of the immune cell infiltration results was set to 0.05, and the vioplot package was used to draw the difference violin diagram and the ggpubr package to draw the boxplot. Differential expression of immune cells including naïve B cells, B cells memory, T cells CD4 memory activated/resting, dendritic cells activated/resting, eosinophils, macrophages M0-2, mast cells activated/resting, monocytes, neutrophils, NK cells activated/resting, T cells CD4 naïve/CD8/follicular helper/gamma delta/regulatory (Treg) by clustering group. Fraction difference, estimate score, immune score, and stromal score were assessed by clustering group as well. Then, GSEA was conducted by GSEA software to explore possible metabolic pathways these lncRNAs were involved in.

LASSO regression was performed by R package glmnet and calculated the minimum criteria. Then, TCGA data were divided into training and testing sets, which were used for survival analysis, survival receiver operator characteristic (ROC) curve drawing, risk score, and risk score related heatmap.

Clinical factors, including age, gender, grade, stage, and risk score in both training and testing sets were computed hazard ratio (HR) of survival by both univariate and multivariate analysis. After that, differential expression and survival analysis by clinical factors were performed using low- and high-risk score group. The model validation of clinical grouping was performed with the survival and survminer packages, and the p-value was used to observe whether the model was suitable for patients with different clinical traits, and the ggpubr package was used to draw the heatmap and boxplot of the clinical traits. Moreover, correlation plots between risk score and immune cell as well as differential expression between low- and high-risk score groups were drawn.

Finding the miRNAs regulated by AC012073.1 from the Starbase (https://starbase.sysu.edu.cn/) website, and selected the top five microRNAs according to agoexpnum from large to small. The downstream mRNA was searched according to miRTarBase (https://mirtarbase.cuhk.edu.cn/~miRTarBase/miRTarBase_2022/php/index.php) and TargetScan (http://www.targetscan.org/vert_71/). According to the corresponding regulatory relationship, Cytoscape software was used to draw the lncRNA-miRNA-mRNA interaction map. For the obtained mRNA targets, prognostic correlation analysis was performed: merged the files of survival data and expression data, used the survival R package to perform Cox analysis, set the significance filter condition p value to 0.00001, and obtained the prognostic mRNA and the fluctuation range of HR values. Using the plot function to draw a forest graph, set high-risk mRNAs to be represented in red, and low-risk mRNAs to be represented in green. The perl script was used to extract the expression data files and phenotype data files of the target mRNA from the gene expression matrix, and the GSEA enrichment analysis was performed with the GSEA software. Using the same method, the downstream target prediction analysis of all prognosis-related lncRNAs was performed, and the downstream regulated microRNAs and mRNAs were obtained, and the ceRNA network interaction map was drawn according to the regulatory relationship. According to the interaction network, mRNA targets co-regulated by two or more lncRNAs were selected as potential co-acting targets.

A total of 10 pairs specimen (10 specimen of Barcelona Clinical Liver Cancer (BCLC) C stage as experimental group, 10 specimens of BCLC A stage as control group). All patients signed informed consent and the study was approved by the Ethical Committee of The First Affiliated Hospital of Guangxi Medical University (2021(KY-E-336)). Specimens were collected after surgery with approximately 100 mg and stored into EP pipe with 1 ml pre-cooled RNAstore Reagent (cat. No. DP408-02; TIANGEN BIOTECH). Then, specimens were stored in -80°C refrigerator.

Cell lines of Huh7 and Hep3B were purchased from Procell Company and cultured in the DMEM (Dulbecco’s Modified Eagle Medium; cat no.8121587; Gibco) culture with 10% FBS (fetal bovine serum; cat no.2146373; VivaCell) and 1% penicillin and streptomycin (cat no.P1400; Solarbio). Total RNA of specimen and cell line were extracted by NucleoZOL reagent (cat no.740404.200; Gene Company Limited). Reverse transcription was used (TransScript^® Uni All-in-One First-Strand cDNA Synthesis SuperMix) for qPCR (One-Step gDNA Removal) (cat no.AU341-02; Transgen). Roche LightCycler 480 machine was employed to perform RT-PCR. FastStart Universal SYBR Green Master (Rox) (cat no.04913914001; Roche) was used for RT-PCR. GAPDH was used as inner control and the relative expressions of lncRNA LOC105377626 were calculated by 2^-△△CT method. Primers of GAPDH and LOC105377626 were as follows: (GAPDH) Forward: GCACCGTCAAGGCTG AGAAC, Reverse: TGGTGAAGACGCCAGTGGA, (LOC105377626) Forward: TCTGTGCTGGACTCCTGCTACTC, Reverse: TGCTTCGGATAGAGGTGCTGAGG. The reaction conditions are set as follows: Hold stage: 50°C for 2 min, 95°C for 2 min. Amplification reaction stage: 95°C for 15 s, 60°C for 1 min, a total of 40 cycles. Melt stage: 95°C for 15 s, 60°C for 1 min, 95°C for 30 s, 60°C for 15 s. The melting curve and amplification curve of each gene in the whole reaction process are automatically generated by the instrument.

A total of 5 x 10⁶ cells were collected and total RNA was performed according to the manufacturer’s instructions (cat no.01030165; Invitrogen, Ambion^® PARIS™). Reverse transcription and RT-PCR assays were performed according to the above procedures. Then, LOC105377626 expressions in both cytoplasmic and nuclear were calculated separately to distinguish cellular location of LOC105377626. Cells were seeded into 6-plate wells and then transfected with antisense oligonucleotide (100 nM, Ribobio Company) and Lipo6000 reagent (cat no.C0526; Beyotime Biotechnology) in 500 ul DEME volume for 24h. After 48h, cells were digested for further experiments.

Transfected cells (2000 cell/pore) were seeded into 6-plate wells for colony formation assay. Twenty days later, cells were fixed by 4% paraformaldehyde (cat no. P1110; Solarbio) for 15 min and dyed by 0.1% crystal (cat no.G1064; Solarbio) violet for 10 min. Then, colonies were counted and compared between groups. Transfected cells (1000 cell/pore) were seeded into 96-plate wells for CCK-8 assay. Cells were cultured with a mixture of 10ul of CCK8 reagent (cat no. MA0218; Meilunbio^®) and 100ul of serum-free DMEM medium for 1h and measured OD values of 450 nm at 0h, 24h, 48h, 72h, 96h, and 120h. Then, proliferation curves were drawn and differences were compared between groups. Transfected cells (3 x 10⁴ cell/pore) were seeded into 96-plate plates for EdU assay (Cell-LightTM EdU Apollo567 In Vitro Kit; cat no.C10310-1; Ribobio Company). Procedures were conducted according to the manufacturer’s instructions.

In migration assay, transfected cells (1x10⁵cell) were seeded into upper room of Transwell plate (8μm pore size, Corning Costar, USA). Lower room was added 700 μl DMEM with 20% FBS and Transwell plate was cultured for 24h in the incubator at 37°C. Transwell plate was fixed by 4% paraformaldehyde and dyed with 0.1% crystal violet. Then, cells were counted in the microscope and compared between groups. In invasion assay, Transwell plates were added 60 μl Matrigel (cat no.354248; Corning Costar, USA) after the Matrigel were diluted 10 times with DMEM. Then, the Transwell plate was placed into the incubator at 37°C for 1h. The rest of the procedure was conducted the same as in the migration assay. Transfected cells (3 x 10⁵ cell/pore) were seeded into 6-plate wells for 24h. Then, scratch of each pore with “+” was conducted and cells were cultured in DMEM with 1% FBS. After 48h, cells were pictured by the microscope and distances and migration distance were measured for comparison.

A P ≤ 0.05 was identified as a significant statistical difference in comparing groups. Non-paired T test or Mann Whitney Test was used to compare the differences between two groups. Differential expression, survival analysis, and interaction network were performed by R (version 4.1.1, https://www.r-project.org/) platform. For Kaplan-Meier curves, p-values and hazard ratio (HR) with 95% confidence interval (CI) were generated by log-rank tests and univariate Cox proportional hazards regression. Cox analysis was used to analyze survival-related risk factors. If the HR value was greater than 1, the factor was a risk factor; if HR<1, the factor was a protective factor; if HR=1, the factor had nothing to do with patient prognosis. The confidence interval of the Cox analysis was set as 95%. Unpaired T test was used to compare the results between groups in the in vitro cell experiment; Mann Whitney Test was used to compare the differences in lncRNA expression in the tissues of patients with stage A and stage C.

A flow diagram was drawn to exhibit the analyzing process of the present study ( Figure 1A ). Then, interaction network between m6a modification-related mRNAs and lncRNAs was established in Figure 1B . Among 23 m6a modification-related mRNAs, 21 of them were enriched in the network, such as IGFBP2/3, YTHDC1/2, METTL3/14/16, and YTHDF1/2/3. After that, univariate survival analysis found 26 prognosis-related lncRNAs, with all of them indicating higher risk of survival (all P<0.001, HR>1, Table 1 ), including AL158166.1, SNHG3, LINC01224, AL442125.2, DUXAP8, NRAV, etc. In addition, differential expression identified 14 lncRNAs, such as SNHG3, AC026256.1, AL031985.3, SNHG4, NRAV, etc., all high in tumor tissues and heatmap visualized their expressions ( Figures 1C, D ). Interestingly, 14 differential lncRNAs were consistently shown in 26 prognosis-related lncRNAs.

Clustering was divided into eight groups labeled2 through 9. According to consensus cumulative distribution function and delta area ( Figures 2A, B ), Cluster 2, Cluster 3, and Cluster 4 drew the survival plots ( Figures 2C–E , all P<0.001) and only two of the groups were picked for further analysis concerning patient numbers in each group. Then, differential expression of CD274 was analyzed and found only in tissue type but not cluster ( Figures 2F, G , P<0.001, P>0.05). Then, expression and correlation heatmaps were established to indicate correlations between clinical factors and PD-L1 with these survival-correlated lncRNAs ( Figures 2H, I ). PD-L1 was correlated with 11 of these lncRNAs ( Figure 2I , all P ≤ 0.05).

Differential analysis was conducted in immune cells by cluster group. Among them, three types of immune cells demonstrated differential expression, with a higher level in naïve B cells and T cells CD4 memory resting and lower level in B cells memory ( Figures 3A–C , 4A , all P<0.05). Differential expression did not show in other types of immune cell ( Figures 3D–V , all P>0.05). Differential score was assessed and only showed a difference in stromal score but did not ESTIMATE score and immune score ( Figures 4B–D , P=0.34, 0.62, 0.018). On the other hand, exploration of possible mechanisms by GSEA showed that these lncRNAs participated in spliceosome, cell cycle, base excision repair, RNA degradation, ribosome, Notch signaling pathway, DNA replication, purine metabolism, etc. ( Figures 4E–L , all P<0.05, false discovery rate<0.25).

LASSO regression was performed and identified the minimum criteria ( Figures 5A, B ). It has been applied for the optimal choice to avoid overfitting as well as to discern the most meaningful model for predicting survival outcome. Then, 370 HCC patients were divided into training set (186 patients) and testing set (184 patients). Survival analysis of two sets was conducted and both indicated significant differences ( Figures 5C, D , P<0.001, P=0.037). Then, survival ROC curves and risk score models, including risk score, survival status, and heatmap, were calculated and visualized, with both sets showing high AUC (training: 0.814, testing: 0.779, Figures 5E–H ).

Risk score was assessed by clinical factors and showed a different risk score in grade, T, and cluster( Figures 6C, D, G , all P<0.05), but did not in age, gender, N, M, and immunesocre ( Figures 6A, B, E, F, H , all P>0.05). Then, heatmap was established to indicate lncRNA expression, clinical factors, risk score, immune score, and cluster ( Figure 6I ). Moreover, both univariate and multivariate analysis were performed and found stage and risk scores steadily influenced clinical outcomes (all P<0.05, HR>1) but not for age, gender, and grade in training and testing sets ( Figures 6J, K ). After that, survival analysis by differential clinical factors was performed and demonstrated a significant difference in age (<65, >65 years old), gender (male, female), G (G1-2, G3-4), T (T1-2, 3-4), stage (I-II, III-IV), M0, and N0 ( Figures 7A–G ) but did not in M1 group. PD-L1 differential expression analysis found that the high-risk group showed higher PD-L1 expression ( Figure 7H ,P=0.043). Correlation analysis demonstrated that T cells regulatory was positively correlated to risk score ( Figure 7I , r=-0.27, P=0.044).

A total of five miRNAs regulated by lncRNA AC012073.1 were obtained by prediction, of which three miRNAs regulated 39 mRNAs, and an interaction network diagram was drawn according to their regulatory relationships ( Figure 8A ). Cox regression was used for prognostic analysis, and a total of 17 prognostic-related mRNAs were screened out, all of which had p-values less than 0.01 and were all high-risk genes ( Figure 8B ). KEGG enrichment analysis showed that these genes were mainly involved in biological activities such as SPLICEOSOME, PYRIMIDINE_METABOLISM, RNA_DEGRADATION, PURINE_METABOLISM and CELL_CYCLE ( Figures 8C, D ). The downstream targets of prognosis-related lncRNAs were predicted, gene targets that were jointly regulated by more than two lncRNAs were screened, and two target genes, MYLIP and HDGF, were obtained ( Supplemental Figure S1 ).

Then, sequence analysis found that AC012073.1 includes an uncharacterized sequence, namely LOC105377626. After removing values that deviate too much, differential expression found that BCLC stage C had higher levels than BCLC stage A ( Figure 9A , P=0.0079). Assay of RNA cytoplasmic and nuclear separation suggested that LOC105377626 was mainly congregated in nucleus compared to cytoplasm (Hep3B: foldchange=19.875, Huh7: foldchange=5.238, Figure 9B ). After that, LOC105377626 was down-regulated by antisense oligonucleotide ( Figure 9C , P ≤ 0.001). We used colony formation, CCK-8, and EdU assays and found that down-regulated LOC105377626 demonstrated reduced colony formation ability ( Figure 9D , foldchange=2.148, P ≤ 0.01) and growth ability ( Figures 9E, F , foldchage=1.545, all P ≤ 0.05). Moreover, Transwell (with and without Matrigel) and scratch assays also consistently suggested that down-regulated LOC105377626 demonstrated reduced invasion ( Figure 10A , P ≤ 0.05) and migration abilities ( Figures 10B, C , P ≤ 0.05).