SMMC-7721 cells were obtained from Cell Bank of Type Culture Collection of Shanghai Institute of Biochemistry & Cell Biology, Chinese Academy of Science. They were cultured in Dulbecco’s Modified Eagle’s Medium (Gibco, California, USA) supplemented with 10% fetal bovine serum (FBS) (Gibco). SMMC-7721 cells were identified by their morphological characteristics which were consistent with the establishment report (12). Cells were not contaminated by mycoplasma or infected with bacteria or fungi. All cells were cultured in a humidified incubator with 5% CO₂ at 37°C. For the sulfatide treatment, cells were incubated at the initial density of 0.5×10⁵ cells/mL and treated with 2 μM galactocerebroside (Gal-Cer) or sulfatide (Sigma, St. Louis, MO, USA).

Microarray profiling was conducted in the laboratory of Aksomics Inc. (Shanghai, China). The microarray was analyzed using the nrStar™ Functional LncRNA PCR chip software, version 1.0 (ArrayStar, Rockville, MD, USA). The hierarchical clustering analysis was carried out using a platform-independent software TBtools (version x64_1_09867) (13).

The TCGA Liver Cancer project (TCGA-LIHC) (N=423) data were downloaded from the UCSC Xena database (https://xenabrowser.net/) (14). Log2(x+0.001) transformation was used to standardize every gene expression profiles and noncoding genes were identified based on their Ensemble gene IDs.

A Cox proportional-hazards regression model was established to analyze the relationship between sulfatide-related lncRNA expression and overall survival (OS) in HCC. The patients were divided into two groups according to the best cutoff value for each sulfatide-related lncRNA, which calculated by the R package maxstat. The OS significance map in HCC was evaluated using the Kaplan-Meier plotter (http://kmplot.com/analysis/) (15).

The least absolute shrinkage and selection operator (LASSO) regression analysis was conducted on the sulfatide-related lncRNAs. The LASSO regression algorithm was used for feature selection with 10-fold cross-validation. The R package glmnet was used for the analysis. For Kaplan-Meier curves, p-values and hazard ratios (HRs) with 95% confidence intervals (CIs) were generated using log-rank tests and univariate Cox proportional -hazards regression. Finally, six sulfatide-related lncRNAs were selected for incorporation into the risk score. The regression coefficient β for multivariate Cox regression model and lncRNA expression were used to construct the risk score formula as follows:

The Tumor Immune Estimation Resource (TIMER) database (http://timer.cistrome.org/) (16) analyzes immune cell infiltration in tumor tissues using high-throughput sequencing (RNA-Seq expression profiling) data (16, 17). The B cell, CD4+ T cells, CD8+ T cells, neutrophil, macrophage and dendritic cells infiltration score of HCC are evaluated by the Timer method of IOBR (version 0.99.9) (18), an R software package, based on the expression profile data of TCGA-LIHC. Spearman’s correlation coefficient between NRSN2-AS1 and immune cell infiltration score in HCC was calculated using corr.test function of R package psych (version 2.1.6) to determine the significantly correlated immune infiltration score.

Total RNA samples were extracted, ragmented to 100bp, and immunoprecipitated using anti-m6A antibody (abcam). Then, eluted RNA and MeRIPed RNA were analyzed using deep sequencing with an Illumina Novaseq™ 6000 platform on the CLOUDSEQ Bio-tech Ltd (Shanghai, China) following the vendor’s recommended protocol (19)l.

Pathway enrichment analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) was performed using the Gene Set Enrichment Analysis (GSEA) database in TCGA-LIHC, which classified the data into high- and low- expression groups based on their NRSN2-AS1 expression (20). Gene sets with |normalized enrichment score (NES)|>1, nominal p value (NOM P) <0.05, and false discovery rate (FDR) q <0.25 were considered to have significant enrichment.

Differences in the expression of sulfatide-related lncRNAs between normal and tumor samples from each tumor were analyzed for significance using unpaired Wilcoxon rank sum and signed rank tests. Survival curves were statistically tested using the log rank test, where p-values and HRs with 95% CIs were represented via Kaplan-Meier plots. The significant correlations between sulfatide-related lncRNAs expression and immune cell infiltration scores in HCC were determined by analyzing Spearman’s correlation coefficients. Pearson correlation analysis was performed between the expression level of NRSN2-AS1_(ENSG00000225377) and gene set expression level of RNA m6A methylation-modifying enzyme. The level of significance was set at P < 0.05. The bioinformatics analysis platform Sangerbox, version 3.0 (http://vip.sangerbox.com/) (21), was used for processing of all the statistical analyses.

The lncRNA profiles of sulfatide-treated HCC cells were compared to those of control cells using ArrayStar lncRNA microarray V2.0. This comparison identified 85 differentially expressed lncRNAs (DE-lncRNAs) based on their Ensemble IDs (|FC|>2, P < 0.05) ( Figures 1A, B . Supplementary Table 1 ). These DE-lncRNAs were further classified by biotype, most of which were lncRNAs as processed transcripts, unclassified processed transcripts, processing/unprocessed pseudogenes, and small amounts of protein coding transcripts or unclassified transcripts ( Figure 1C ).

The TCGA-LIHC dataset was used to detect the expression of these 85 sulfatide-related lncRNAs in HCC tissues, showing that 24 of them were highly expressed in HCC compared to normal liver tissues ( Figure 2A ). However, the expression of three sulfatide-related lncRNAs, AP002841.2 (ENST00000504610), RP11-733O18.1 (ENST00000422914) and RP5-885L7.10 (ENST00000412500) were lower in HCC tissues ( Figure 2B ).

Among the above mentioned 27 sulfatide-related lncRNAs, six lncRNAs were ultimately identified to be related to prognosis ( Figure 3A ), including RP11-122M14.1 (ENST00000415202), RP11-280O1.2 (also known as LRRC52-AS1; ENST00000438275), AC079354.5 (ENST00000447111), AC005037.3 (ENST00000413848), AC108488.3 (also known as RNASEH1-AS1; ENST00000438436) and RP5-1103G7.4 (also known as NRSN2-AS1; ENST00000442637). Kaplan–Meier survival analysis was utilized to evaluate the significance of lncRNA expression in patient prognosis ( Figure 3B ). High levels of these sulfatide-related lncRNAs were all correlated with poor prognosis in patients with HCC ( Figure 3 ).

Based on the expression of six sulfatide-related lncRNAs and multivariate Cox regression coefficient, the prognosis risk score for sulfatide-related lncRNAs was calculated using the following formula: riskscore (lambda.min=0.0029) = (2.0727) × LRRC52-AS1 + (0.3691) × RNASEH1-AS1 + (0.2646) × NRSN2-AS1 ( Figures 4A, B ). Subsequently, an X-tile diagram was used to generate the optimal cutoff point for the risk score. The TCGA-LIHC patients were divided into high- and low-risk groups based on this cutoff risk score value. A prognostic curve and a scatter plot were used to indicate the risk score and survival status of each HCC patient ( Figures 4C, D ). In addition, the heat map of the expression profiles for candidate lncRNAs demonstrated that they were all highly expressed in the high-risk group ( Figure 4E ). Kaplan-Meier analysis validated that the TCGA-LIHC patients in the high-risk group showed a significantly worse survival than those in the low-risk group at the 10-year time point ( Figure 4F ). Furthermore, the time-dependent receiver operating characteristic (ROC) analyses showed that the area under curve (AUC) for the risk score model was 0.671 at the one-year time point, 0.621 at the three-year time point, and 0.629 at the five-year time point ( Figure 4G ). Taken together, these findings represented the three sulfatide-related lncRNAs as the prognostic signature for HCC patients.

Our previous study reported that sulfatide does not only affect the binding of METTL3 to METTL14 and WTAP by acetylating the METTL3 protein (19), but also inhibits the YTHDF2 expression in HCC cells (22). Next, we investigated whether RNA m6A methylation modification affected sulfatide-related lncRNA expression. The MeRIP-seq experiments were performed, in order to clarify the role of RNA m6A methylation modification. Their results showed that the abundance of m6A in NRSN2-AS1, one of sulfatide-related lncRNAs, was significantly increased in sulfatide-treated HCC cells. This suggested that m6A modification was related to the regulation of NRSN2-AS1 expression ( Figure 5A ). Furthermore, the relationship between NRSN2-AS1 expression level and a series of m6A-binding proteins (23) was analyzed in TCGA-LIHC samples. The results showed that the expression levels of NRSN2-AS1 were positively correlated with the expression of m6A writer and reader signatures ( Figure 5B ). In summary, the expression of NRSN2-AS1 in HCC was related to the changes in RNA m6A methylation induced by sulfatide.

An increasing number of studies have demonstrated that sulfatide is involved in tumor immunity, where the HIF-1-galactose-3-o-mercaptotransferase 1-sulfide axis enhanced immune escape of renal clear cell carcinoma by increasing tumor cell-platelet binding (24). In addition, a subpopulation of type II Natural killer T cells (NKT cells) characterized by their response to autoglycolipid sulfides was shown to induce a major immunomodulatory mechanism that controls inflammation in anticancer immunity (25). To further examine whether the sulfatide-related lncRNA NRSN2-AS1 can act as an immune indicator, a correlation analysis of NRSN2-AS1 expression with immune infiltration was performed. TIMER data showed that high NRSN2-AS1 expression was significantly associated with six types of immune cells (B cells, CD4+ T cells, CD8+ T cells, macrophages, neutrophils and dendritic cells) in HCC ( Figure 6 ). This result pointed out that NRSN2-AS1 may serve as an indicator in tumor immune microenvironment (TIME) characterization in HCC.

To investigate the biological functions and pathways associated with the sulfatide-related lncRNA NRSN2-AS1, the TCGA-LIHC samples were divided into high- and low-expression groups based on their NRSN2-AS1 expression. GSEA was used to evaluate the enrichment of KEGG pathways. The pathways associated with high NRSN2-AS1 expression were enriched in the Cell Cycle pathway ( Figure 7A ). The pathways associated with low expression of NRSN2-AS1 were enriched in peroxisome and perxisome proliferator-activated receptor (PPAR) signaling pathways related to immune response (26, 27), as well as a variety of amino acid (tryptophan, arginine, proline, glycine, serine, threonine, tyrosine, and histidine) and lipid (fatty acid and linoleic acid) metabolic pathways ( Figures 7B-D ).