Gene expression data and clinicopathological characteristics of HCCs used in this study were generated by the Cancer Genome Atlas (TCGA) and are available at https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga. Downloaded data included FPKM (fragments per kilobase of transcript per million) reads-based gene expression data and the raw read count values. The R package “TCGAbiolinks” was used for downloading (22). After screening for data quality, 374 HCC samples were retained in this study. Of the 374 cases, one lacked prognostic information, so 373 were used for model construction and survival analyses (Figure 1). In addition, the tumor immune dysfunction and exclusion (TIDE) algorithm, as described in a previous study (15), was used to determine both the T cell exclusion (TCE) level and the T cell dysfunction level using the FPKM expression matrix.

To explore what chemokines/cytokines and immune checkpoint ligands mediate the communications between HCC cells and CD8+ T cells, we analyzed a single-cell RNA-sequencing (scRNS-seq) data of HCC (GSE146115) available in the Gene Expression Omnibus (GEO) (23). After performing imputation on the dropouts by the “scImpute” algorithm (24), the R package “Seurat” (25) was used for dimensional reduction, clustering analysis, and cell type annotation. Finally, the CellChat Explorer (26) program was used to infer the biologically significant interactions between chemokines, cytokines, and immune checkpoint (ICP) ligands and their receptors in the interactions between HCC cells and CD8+ T cells.

To identify lncRNAs that are differentially expressed between HCC and adjacent morphologically normal liver tissues, we used the R packages “edgeR” (27, 28) and “limma” (29) to analyze the 374 HCC tissues and 50 adjacent normal liver tissues. The thresholds of P ≤ 0.05 and |log2FC| (FC: Fold change) > 0.5 were used. The GENCODE database (30) was used to identify lncRNAs.

After HCC tissues were divided into TCE-higher (n = 187) and TCE-lower (n = 187) groups by the median TCE level following the TIDE analysis, the edgeR-limma procedure was also used to identify differentially expressed lncRNAs between the two TCE groups.

LncRNAs differential expression between HCC and normal tissues and between TCE-higher and TCE-lower groups were identified as HCC- and TCE-associated lncRNAs.

The 373 HCCs with prognostic information were randomly assigned to the training cohort (n = 187) and validation cohort (n = 186) at a 1:1 ratio using the R package “caret”. Univariate Cox regression analysis was performed to assess the association of each differentially expressed and TCE-associated lncRNA with the overall survival (OS) in the training cohort. LncRNAs significantly correlated with OS (P ≤ 0.001) were subjected to the LASSO-COX regression analysis (31) to develop the prognostic model (i.e., 11LNCPS). Based on the model, a risk score (RS) for OS was built based on a linear combination of the regression coefficient derived from the multivariate Cox regression model and the expression level of the optimized lncRNAs.

The 11LNCPS score (risk score) was computed as follows: Risk score = ∑ i = 1 N ( C i × F a c t o r i ) , where N represents the number of prognostic factors, Factor_i represents the expression of lncRNAs, and C_i represents the regression coefficient of the multivariate Cox regression model (1, 32, 33).

HCCs in the training cohort were then divided into two groups using the median, one with higher 11LNCPS scores and the other with lower scores, for model evaluation. Kaplan-Meier analysis was used for overall survival, with the log-rank test to evaluate statistical significance.

The time-dependent receiver‐operating characteristic (ROC) analysis was performed to evaluate the accuracy of the 11LNCPS classifier in survival prediction at 1, 2, and 3 years in the training cohort. The calibration curves, Harrell’s concordance index (C-index) curves, and ROC Areas under the curves (AUCs) were calculated for model evaluation. Calibration curves were calculated by calibrating the function implemented in the “rms” package to assess the predictive ability. The 11LNCPS was also validated in the validation and entire cohorts. In addition, the 11LNCPS model was compared with two previously reported effective HCC prognostic models for ROC, C-index, and prediction error curves in the validation cohort using the packages of “MASS”, “timeROC”, “survival”, and “survminer”. One was the 8-gene model containing H2AFX, SQSTM1, ITM2A, PFKP, TPD52L1, ACSL4, STRN3, and CPEB3 (34); and the other was the 4-gene model containing CENPA, SPP1, MAGEB6, and HOXD9 (35).

Considering that TCE is primarily related to the immune escape (15), we applied the xCell computational method (36–38) to estimate the enrichment scores (xCell scores) of different immune cell types in HCCs with higher and lower 11LNCPS scores.

A total of 374 HCC samples with normalized gene expression FPKM data and standard annotation were used to analyze the distribution of 34 types of immune cells using the xCell pipeline. The 34 types of immune cells included CD4+ naive T-cells, CD4+ T-cells, CD4+ memory T-cells, CD4+ Tcm (central memory T cell), CD4+ Tem (effective memory T cell), CD8+ naive T-cells, CD8+ T-cells, CD8+ Tcm, CD8+ Tcm, Treg cells, gamma delta T cells (Tgd cells), Th1 cells, Th2 cells, natural killer T cell (NKT), natural killer cell (NK), pro B-cells (B cell progenitors), naive B-cells, B-cells, memory B-cells, class-switched memory B-cells, plasma cells, monocytes, macrophages, macrophages M1, macrophages M2, dendritic cell (DC), activated dendritic cell (aDC), conventional dendritic cell (cDC), plasmacytoid dendritic cell (pDC), immature dendritic cell (iDC), neutrophils, eosinophils, mast cells, and basophils.

The Kaplan-Meier (K-M) analysis was performed to assess whether the infiltrating status of different immune cell types affects patient survival in HCC. Survival outcomes were calculated and visualized using the R packages “survival” and “survminer”. The correlation between an xCell score and an immune cell type was analyzed with P ≤ 0.05 and |log2 (FC)| > 0.25 following the procedure described previously (36), and the outcome was visualized using the R packages “pheatmap”, “EnhancedVolcano” and “ggpubr”.

TIDE algorithm was then applied to evaluate the association of 11LNCPS scores with TCE and T cell dysfunction. HCCs were divided into higher and lower 11LNCPS scores using the median, and the TIDE algorithm (15) was then applied to each group. The outcome was visualized using “ggpubr”.

To determine whether the 11LNCPS score is associated with therapeutic responses to ICIs, the 373 HCCs were divided into two groups, one with higher and one with lower 11LNCPS scores, using the median. The Subclass Mapping (SubMap) algorithm (39) was then applied to measure the correspondence between the two 11LNCPS groups and groups of malignancies with and without responses to anti‐CTAL‐4, anti‐PD‐1, and anti‐PD‐L1 therapies from previous studies (40, 41). The outcome was visualized using the R packages of “pheatmap” and “ggpubr”.

To determine the critical biological pathways and characteristics of HCC determined by 11LNCPS score, the GO, KEGG pathway analysis, and GSEA (42, 43) were applied to HCCs with higher and lower 11LNCPS scores using the R packages “GSEABase”, “clusterProfiler” (44), “enrichplot”, and “org.Hs.eg.db”. Briefly, the edgeR-limma procedure was used to find differentially expressed genes (DEGs) between HCCs with higher-RS and those with lower RS in TCGA. DEGs with thresholds of P ≤ 0.05 and |log2FC| (FC: Fold change) > 0.5, were subjected to GO and KEGG analysis. Go analysis included BP (biological process), CC (cellular component), and MF (molecular function). For GSEA, all DEGs were subjected to the “GSVA” package (45) after ranking from high to low based on their FC values. A P value smaller than 0.05 was considered significant in the GSEA. The hallmark gene set “h.all.v7.1.symbols.gmt” was downloaded from https://www.gsea-msigdb.org/and subjected to “GSVA” in a similar fashion.

For each of the 11 lncRNAs in the 11LNCPS, a series of analyses were performed to identify the core one. The Kaplan-Meier survival analysis was performed to assess whether a lncRNA’s expression level is associated with the overall survival (OS) in the 373 HCC patients.

Correlation of each lncRNA’s expression level with T cell exclusion (TCE) was ranked based on the Spearman correlation coefficient value, with those greater than 0.2 with P ≤ 0.05 considered significant. Infiltration levels of prognosis-associated immune cells were also compared between HCCs with higher and lower 11LNCPS scores and HCCs with higher and lower expression levels of the 11LNCPS lncRNAs. Those with a lower level. LncRNAs whose higher expression levels significantly correlated with worse patient OS, whose Spearman correlation coefficient values were greater than 0.2 (P ≤ 0.05), and that affected infiltrations of more immune cell types were considered crucial members of the 11LNCPS, including LINC01134 and AC116025.2. The outcome was visualized by the R packages “pheatmap”, “ggpubr”, “corrplot” and “ggplot2”.

The relationship between LINC01134 and AC116025.2 expression and TCE or T cell dysfunction was tested using the TIDE algorithm, and the outcome was visualized using the “ggpubr” R package.

To explore the biological functions of LINC01134 and AC116025.2 in HCC, we performed GO, KEGG, and GSEA analysis in HCCs as described in the previous enrichment analysis.

Normal liver cells QSG-7701 and LO2 were kindly provided by Dr. Liang Yang of the Southern University of Science and Technology. HCC cell lines HepG2 and Huh-7 were purchased from the BeNa Culture Collection (Beijing, China). The Jurkat cell line was kindly provided by Dr. Lili Ren of Shenzhen People’s Hospital. The DMEM medium (Gibco, USA) supplemented with antibiotics (Biological Industries, Israel) and 10% FBS (Gibco) was used for liver cell culture. The RPMI 1640 medium (Gibco) supplemented with antibiotics and 10% FBS were used for Jurkat cells. All cells were cultured at 37 °C in a humidified atmosphere containing 5% CO₂.

Both the negative control siRNAs (si-NC) and the LINC01134 siRNAs were provided by GenePharma (Shanghai, China). Sequences of siRNAs against LINC01134 were 5′-GACAGGTTTGAGCTAGAAAC-3′ (si-LINC01134-1) and 5′-GCAAAUGCACAGCGAGGAAAG-3′ (si-LINC01134-5). At confluency of 30–50%, HepG2 or Huh-7 cells were transfected with siRNAs using the Lipofectamine RNAiMAX reagent (Invitrogen, USA). After 48 hours, transfected cells were split into two portions. One was used for RNA isolation and gene expression analysis, and the other was grown in a 6-well plate for 48 hours to collect the conditioned medium (CM). Each experiment was repeated twice unless otherwise stated.

Total RNA was extracted from cultured cells using the Eastep Super Total RNA Extraction Kit (Promega, USA) and reverse transcribed into cDNA using the HiScript III All-in-one RT SuperMix Perfect for qPCR Kit (Vazyme, China). PCR was performed with the KT SYBR qPCR Mix (Ktsm-life, China) using the qTOWER 3.0 PCR system (Jena Industries, Germany). Primers and their sequences are as follow: LINC01134, 5′-ATGAACAGCAAATGCACAGCG-3′ (forward) and 5′- ATAGGTCTTGGCTGGTTCTCG-3′ (reverse); AC116025.2, 5′-TGGAGCAGAAAGAGCTGTCTCAAG-3′ (forward) and 5′-TGTCAGGAAACTGTGTGGACG-3′ (reverse); CXCL1, 5′-CTGGCTTAGAACAAAGGGGCT-3′ (forward) and 5′-TAAAGGTAGCCCTTGTTTCCCC-3′ (reverse); CXCL2, 5′-CCCATGGTTAAGAAAATCATCG-3′ (forward) and 5′-CTTCAGGAACAGCCACCAAT-3′ (reverse); CXCL3, 5′-CGCCCAAACCGAAGTCATAG-3′ (forward) and 5′-ACCTTGCCTTCTTTGTCTTTGTTGGA-3′ (reverse); and β-actin, 5′-TCCCTGGAGAAGAGCTACGA-3′ (forward) and 5′-GCTCCCCTTGTTCAGTATCTTTT-3′ (reverse). In the PCR, β-actin served as the endogenous control. The relative expression of genes was calculated using the −2ΔΔCt method.

T cell migration was analyzed using the transwell assay as previously described (46–48). Briefly, Jurkat cells (10⁶ cells/ml) were washed with PBS and serum-starved for 3 hours, 10⁵ cells in 0.1 ml were then seeded onto an 8.0-μm pore size insert (Corning, USA), and 400 μl complete medium or CM were then added to the lower chambers of a 24-well plate (Corning). After incubation at 37°C for 16 hours in an incubator, migrated cells in the lower chambers were collected and counted using an automated cell counter (Invitrogen). The numbers of migrated cells in different groups were normalized by the number of cells from the complete medium group.

The R software (version 4.1.1) was used for all statistical analyses and plot drawings except as specifically stated. Patients were randomly grouped using the “caret” R package. The univariate and multivariate Cox proportional hazards regression analyses were performed using the “survival” package. Kaplan-Meier analysis was used for overall survival, with the log-rank test to evaluate statistical significance. Statistical differences between the two groups were assessed using the Wilcoxon test. The grouping basis (the cutoff point) was the median value of each corresponding index.One-way ANOVA with Bonferroni’s multiple-comparisons test was performed for qPCR and T cell migration analysis using the GraphPad Prism (GraphPad Prism 8). P < 0.05 was considered statistically significant unless otherwise stated.

The workflow of the entire study is summarized in Figure 1. In total, the TCGA database contained 374 HCC cases with gene expression data. All the 374 cases were used for the identification of differentially expressed and TCE-associated lncRNAs. One of the 374 cases lacked prognostic information and thus was excluded for model construction and survival analysis. Using the 374 HCCs and 50 cases of noncancerous liver tissues with expression profiling and other information, two groups of differentially expressed genes (DEGs) were identified from a total of 56493 human genes, including lncRNA, other noncoding RNA, and protein coding genes. One group contained 8191 genes that were differentially expressed between HCCs and normal liver tissues, with 6438 upregulated and 1753 downregulated in HCC (Figure 1 and Figure S1A). Among these 8191 DEGs, 2752 were lncRNAs (Figure 1). The other group contained 4127 TCE-associated genes that were differentially expressed between HCCs with higher TCE scores and those with lower TCE scores, including 2914 upregulated and 1213 downregulated in the TCE-higher group (Figure 1 and Figure S1B). Of the 4127 TCE-associated genes, 702 were lncRNAs (Figure 1). In total, 336 lncRNAs were both differentially expressed and TCE-associated in HCC (Figure 2A and Table S1).

Of the 374 HCC cases, one lacked prognostic information and thus was excluded for model construction and survival analysis. The 373 HCCs with survival data were divided into the training (n = 187) and validation (n = 186) cohorts. Each differentially expressed TCE-associated lncRNA in the training cohort was subjected to the univariate Cox regression analysis to evaluate its association with patients’ overall survival (OS). Fifty-four lncRNAs were significantly associated with prognosis (P < 0.001) (Table S1). The LASSO-Cox regression analysis was then performed, in which tenfold cross-validation was applied to overcome overfitting with an optimal λ value of 0.028393 selected (Figure 2B). A combination of 11 lncRNAs had non-zero LASSO coefficients and thus was the most robust prognostic value (Figure 2C). This combination of 11 lncRNAs was named 11 lncRNA prognostic signature (11LNCPS). The 11 lncRNAs included LINC01134, C2orf27A, LINC00501, AC104066.3, AC034229.4, CASC8, FAM225B, AL451069.3, AL161669.3, AC116025.2 and LINC00632.

To determine the 11LNCPS score, the Cox multivariate regression analysis was used to evaluate each of the 11 lncRNA’s contribution to the 11LNCPS (Table 1), which resulted in the following formular for calculating the risk score (i.e., 11LNCPS score) in an HCC: 11LNCPS score = 0.214579 × expression of LINC01134 + 0.019508 × expression of C2orf27A + 1.045738 × expression of LINC00501 + 1.244493 × expression of AC104066.3 + 0.140677 × expression of AC034229.4 + 0.302498 × expression of CASC8 + 7.231505 × expression of FAM225B + 0.089521 × expression of AL451069.3 + 0.226717 × expression of AL161669.3 + 0.224378 × expression of AC116025.2 + 0.371662 × expression of LINC00632.

To test the validity and effectiveness of the 11LNCPS in HCC, we calculated the 11LNCPS risk score for each case in the training, validation, and entire cohorts; divided HCCs in each cohort into the higher- and lower-risk groups using the median 11LNCPS score; and performed a series of analyses (Figures 2D–F and Figures S2A–D). The Kaplan-Meier analysis demonstrated that the OS rate was better in patients with lower 11LNCPS scores than those with higher scores in each cohort (P ≤ 0.05, Figure 2D).

The area under ROC curve (AUC) for 1, 2, and 3 years reached 0.67, 0.7, and 0.74, respectively, in the training cohort; 0.71, 0.69, and 0.64, respectively, in the validation cohort; and 0.68, 0.7, and 0.69, respectively, in the entire cohort (Figure 2E). These AUC curves indicate a reasonable discrimination power of the 11LNCPS in HCC. Additionally, the 11LNCPS’s C-index was greater than 0.60 for 1, 2, and 3 years in each cohort, showing an excellent predictive accuracy of the 11LNCPS (Figure S2A). Furthermore, the calibration curve demonstrated good consistency for 1, 2, and 3 years in each cohort (Figure S2B).

We also compared our 11LNCPS model with two reported models, i.e., the 8-gene model (34) and the 4-gene model (35) in the validation cohort. For each of the 3 time points (1, 2, and 3 years), 11LNCPS showed a higher AUC value (Figure 2F) and a higher C-index (Figure S2C). Each model’s predicted error line overlapped well with the reference line (Figure S2D), demonstrating a lower predicted error rate for each of the 3 models.

We applied the xCell algorithm to the RNA-seq datasets of the 374 HCCs to determine the infiltration levels of 34 types of immune cells (Figure S3). The correlation between an immune cell infiltration and patient OS was evaluated using the Kaplan-Meier analysis (Figure 3A and Figure S4). Altered infiltrations of 7 types of immune cells were significantly associated with OS (Figure 3). Increased infiltrations of CD8+ naive, CD8+ Tcm, CD8+ T, and pDC cells were associated with better OS, while increased infiltrations of Th1, Th2, and pro B cells were associated with a worse OS in HCC (Figure 3A).

To evaluate the relationship between the 11LNCPS and immune responses to HCC, we divided all HCCs into higher and lower 11LNCPS scores using the median and compared the distribution of different immune cell types between the two groups (Figure 3B and Figure S5). HCCs with higher 11LNCPS scores had decreased infiltrations of CD8+ Tcm, macrophages, macrophages M2, aDCs, and cDCs immune cells and increased infiltrations of Th1, Th2, pro B, B, and basophils cells (Figure 3B). In the Kaplan-Meier analysis, alterations in 4 of the 10 immune cell types were significantly associatged with OS (Figure 3A). The 4 alterations included decreased filtration of CD8+ Tcm cells and increased filtrations of Th1, Th2, and pro B cells (Figures 3A, B).

TIDE is a computer program that models the induction of T cell dysfunction in tumors with higher infiltration of cytotoxic T cells and the prevention of T cell infiltration in tumors with lower levels of such cells (15). To further explore the impact of 11LNCPS lncRNAs on immune responses in HCC, we compared HCCs with higher and lower 11LNCPS scores for TCE and T cell dysfunction levels which were analyzed using the TIDE program. HCCs with higher 11LNCPS scores had significantly higher TCE scores and lower T cell dysfunction levels than those with lower 11LNCPS scores (P ≤ 0.05) (Figures 3C, D).

Using the SubMap analysis, we compared HCCs with higher and lower 11LNCPS scores to malignancies with and without responses to immunotherapies from previous studies (41, 49). HCCs with higher 11LNCPS scores were significantly associated with malignancies that respond to a PDL1 inhibitor (Figure 3E, P < 0.05).

Differentially expressed genes were identified in HCCs with higher and lower 11LNCPS scores (Figure S6). Such genes were analyzed to evaluate the effect of 11LNCPS lncRNAs on different biological processes and signaling pathways using the GO, KEGG pathway, and GSEA analyses (Figure 4).

Many biological processes identified in the GO analysis are involved in the cell cycle and DNA replication. These processes included organelle fission, nuclear division, DNA-dependent DNA replication, cell cycle checkpoint, chromosomal region, DNA replication preinitiation complex, single-stranded DNA helicase activity, ATPase activity, and DNA replication origin binding.

In the KEGG pathway analysis, the top-ranked pathways were involved in cell cycle and ligand-receptor interactions, including cytokine and cytokine receptor-related signaling and the viral proteins’ interactions with cytokines and cytokine receptors (Figure 4B). They also included metabolism-associated pathways such as retinol, drug, and xenobiotics (Figure 4B).

The GSEA analysis resulted in similar findings (Figure 4C). Specifically, signaling pathways related to cell cycle and DNA replication were significantly enriched in HCCs with higher 11LNCPS scores, including G2M checkpoint, E2F targets, cell cycle, and UV response containing DNA replication genes. Signaling pathways related to metabolism, immune function, and cell death were significantly suppressed in HCCs with higher 11LNCPS scores, including fatty acid metabolism, bile acid metabolism, IL6-JAK-STAT3 signaling, IFNα response, and apoptosis.

Therefore, the 11LNCPS appears to affect cell cycle signaling pathways, DNA replication, immune function, and cell death.

To rank the 11LNCPS’s 11 lncRNAs for their contributions to the signature, we analyzed them for the association of expression change with OS and immune responses in HCC. In the Kaplan-Meier analysis, the increased expression in 5 of the 11 lncRNAs was significantly associated with worse OS, including LINC01134, AC104066.3, AC034229.4, AC116025.2, and LINC00632 (Figure 5A and Figure S7).

Based on the TCE scores revealed by the Spearman analysis, increased expression in 8 of the 11 lncRNAs was significantly associated with TCE (P < 0.05). These 8 lncRNAs and their Spearman coefficient values were C2orf27A, 0.41; LINC01134, 0.33; AC104066.3, 0.33; LINC00632, 0.31; AC034229.4, 0.29; AC116025.2, 0.26; FAM225B, 0.24; and LINC00501, 0.13, respectively (Figure 5B).

We further evaluate their effects on immune cell infiltration for the 5 lncRNAs whose expression increase was significantly associated with a worse OS.

While expression change in LINC01134 or AC116025.2 significantly affected the infiltrations of 5 immune cell types, expression change in other 11LNCPS lncRNAs altered 3 or 4 types (Figure 5C and Figure S8). Specifically, LINC01134 and AC116025.2 upregulation was significantly associated with increased infiltrations of Th1, Th2, and pro B immune cells but decreased infiltrations of CD8+ naive T and CD8+ Tcm cells (Figure 5C and Figures S8A, D). For the other 3 11LNCPS lncRNAs associated with OS, AC034229.4 upregulation was associated with increased infiltrations of Th1, Th2, and pro B cells (Figure S8C); and higher levels of AC104066.3 and LINC00632 were associated with increased infiltrations of Th2 and pro B cells and decreasing infiltration of CD8+ naive T cells (Figures S8B, E).

Additionally, HCCs with higher LINC01134 or AC116025.2 expression had higher TCE scores and reduced T cell dysfunction levels (Figures 5D, E).

Similar to the analyses of 11LNCPS for its potential impact on biological processes and signaling pathways, we divided HCCs with higher and lower expression levels of LINC01134 or AC116025.2, identified differentially expressed genes, and performed GO, KEGG pathway, and GSEA analyses (Figure S9).

The most enriched processes for LINC01134 in the GO enrichment analysis included cell chemotaxis and chemokine response related biological processes, chromosome related molecular function, receptor-ligand activity, and chemokine binding cellular component (Figure 6A, left).

In the KEGG pathway analysis, LINC01134 upregulation was significantly associated with diverse immune-related signaling pathways, including chemokines/cytokines and their receptors and T and B cell receptors. Some cancer-associated pathways were identified, including PI3K-Akt, Rap1, cell cycle, glioma, and p53 (Figure 6B, left).

In the GSEA analysis, LINC01134 upregulation was associated with the active cell cycle (e.g., E2F targets and G2M checkpoint). It was also associated with cancer pathways (e.g., MYC targets) (Figure 6C, left). On the other hand, LINC01134 upregulation was inversely related to pathways of immune (IFNγ response, IFNα response, IL6-JAK-STAT3, IL2-STAT5), metabolism (bile acid metabolism, fatty acid metabolism), apoptosis, epithelial-mesenchymal transition, UV response, TGFβ, hypoxia, and p53 (Figure 6C, left).

These results support the role of LINC01134 in the cell cycle, cell death, immunity, chemokine expression, and chemotaxis in HCC.

In the GO analysis, AC116025.2-associated genes were primarily enriched in cell division, catabolic metabolism, chromosome and transporter complex, receptor and channel activities, and oxidoreductase activity (Figure 6A, right).

In the KEGG analysis, AC116025.2-associated genes were enriched for pathways in the cell cycle, DNA replication, metabolism, and apoptosis (Figure 6B, right). Multiple metabolic pathways were enriched, including carbon metabolism, retinol metabolism, fatty acid metabolism, tryptophan metabolism, and propanoate metabolism (Figure 6B, right).

In the GSEA, AC116025.2 upregulation was associated with cell cycle activities and cancer-related pathways such as E2F targets, G2M checkpoint, MYC targets, MTORC1 signaling, and mitotic spindle (Figure 6C, right). On the other hand, AC116025.2 upregulation was associated with reduced activities of signaling pathways related to immune, metabolism, and cell death, including IFNγ response, IFNα response, inflammatory response, bile acid metabolism, fatty acid metabolism, apoptosis, UV response, epithelial-mesenchymal transition, KRAS signaling, TGFβ signaling, and estrogen response (Figure 6C, right).

Immune responses often involve cytokines, chemokines, and their receptors. Therefore, we investigated whether expression changes in LINC01134 and AC116025.2 are associated with chemokines, cytokines, and ICP ligands in HCC. In the scRNA-seq data, CD8+ cells could be annotated (Figure S10A). We thus identified the chemokines, cytokines, and ICP ligands synthesized by HCC cells and could mediate CD8+ T cells’ recruitment using the CellChat algorithm (26).

In total, 22 cytokines and chemokines were identified, including CXCL12, CCL5, CXCL16, CCL16, CXCL10, CCL20, IL7, CCL15, CXCL2, IL15, CCL3, CCL4, CXCL8, CXCL9, CXCL11, CXCL1, CCL28, CCL2, CXCL13, CXCL3, CXCL6, and CCL22 (Figure 7A, left). We also identified 26 ICP ligands that could bind to their ICPs, including HLA-A, HLA-B, HLA-C, HLA-E, CD70, PVR, HLA-F, LGALS9, CEACAM1, HLA-DRA, ICOSLG, HLA-DMA, HLA-DPB1, HLA-DOA, HLA-DRB1, CD86, TNFSF15, HLA-DQB1, HLA-DPA1, HLA-DMB, TNFSF4, HLA-DQA1, CD48, HLA-DOB, RAET1E, and RAET1G (Figure 7A, right). Using the Spearman correlation analysis, we found that LINC01134 upregulation in HCC was negatively correlated with the following genes (R_S > 0, P ≤ 0.05): CXCL1, CXCL2, CXCL3, HLA-C, and HLA-E and was positively correlated with LGALS9 (R_S < 0, P ≤ 0.05) (Figure 7B). For AC116025.2, its upregulation was positively associated with CXCL1, CXCL8, CXCL20, and TNFSF4 (R_S > 0, P ≤ 0.05) (Figure 7C).

To test the impact of LINC01134 and AC116025.2 on HCC, we measured their expression in two HCC cell lines using qRT-PCR and found that both LINC01134 and AC116025.2 were significantly upregulated in HepG2 and Huh-7 HCC cell lines compared to normal liver cell lines QSG-7701 and LO2 (Figure 8A). We also knocked down LINC01134 expression in the two HCC cell lines and measured the expression of three cytokines whose expression correlated with LINC01134 in HCC samples. LINC01134 knockdown significantly increased the expression of CXCL2 and CXCL3 (Figure 8B). Consistent with the upregulation of CXCL2 and CXCL3 by LINC01134 knockdown, conditioned medium from HCC cells with LINC01134 knockdown significantly increased the migration of Jurkat T cells (Figure 8C). These findings support the role of LINC01134 in HCC.