Immune‐related lncRNA signature delineates an immune‐excluded subtype of liver cancer with unfavorable clinical outcomes

Abstract Background Long non‐coding RNAs (lncRNAs) play crucial roles in immune regulation and, therefore, may be closely related to the tumor microenvironment (TME). However, there are few studies regarding the relationship between the lncRNAs and the TME in liver cancer. Methods Firstly, we constructed a lncRNA signature based on the top 10 immune‐inversely related lncRNAs obtained from the ImmLnc database and performed disease‐free survival (DFS) and overall survival (OS) analyses for the patients included in the Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA‐LIHC) stratified by the lncRNA signature. Then, we explored the relationship between the lncRNA signature with distinct mutation profiles and the tumor microenvironment (TME). Results The lncRNA signature was successfully constructed and verified by survival analysis. The high lncRNA signature was correlated with a decreased DFS and OS in liver cancer and other two gastrointestinal cancers. The mutation profiles showed that the Lnc_high group had a higher number of mutations on many genes, mostly enriched in p53 and WNT pathways. The TME results showed that the Lnc_high group had the highest proportion (51%) of lymphocyte depletion‐characterized immune subtype, and a higher expression of immune checkpoint molecules such as LAG3, PD‐L1, CTLA4. On the contrary, in the Lnc_low group, infiltrating immune‐cell proportions were significantly higher, and a significant enhancement of four axes of the cancer immunity cycle immunogram was observed in this group. Conclusions The lncRNA signature we constructed identified an immune‐excluded subtype of liver cancer with unfavorable clinic outcomes, which could be tested as a biomarker for immunotherapy in the future.


| INTRODUC TI ON
Liver cancer is one of the most common cancers worldwide with its incidence ranking fifth and mortality ranking fourth among all malignancies in both sexes. 1 For patients with early-stage liver cancer, radical treatment of surgery is the first choice, but most liver cancer patients are diagnosed at an advanced stage. 2 Considering the limitations and poor outcomes of existing treatments, immunotherapy is a promising therapeutic approach; however, only a small percentage of patients benefit from this treatment. 3 Therefore, continued research into molecular biomarkers and novel therapies is critical for predicting prognosis and determining personalized treatment.
The growing data from research on tumor microenvironment (TME) indicates that tumor immune status plays a critical role in cancer progression and prognosis. 4 Multiple TME markers, such as tumor-infiltrating immune cells and immune-related gene expression profiles, have provided extensive evidence in diverse cancer types that these markers are useful, not only in cancer prognosis but also in cancer immunotherapy. 5 Long non-coding RNAs (ln-cRNAs), a type of non-coding RNAs that are greater than 200 nucleotides in length, are critical regulators of gene expression and play a crucial role in immune regulation. For instance, the lncRNA KCNQ1 can play an immunosuppressive role, thereby promoting the immune escape of liver cancer. 6 Recent studies have found that lncRNAs are correlated with immune-cell infiltration in different cancers, including liver cancer. 7 Furthermore, several studies have reported that lncRNA profiles can be used to predict the prognosis of liver cancer. 8 However, there are currently few studies regarding the relationship between the lncRNAs and the TME in liver cancer.
In this study, we analyzed the Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) cohort data and developed a novel immune-related lncRNA signature based on the lncRNAs that significantly affect immune pathway activity. We then investigated the prognostic value of this signature in the cohort and evaluated the relationship of the lncRNA signature with mutation profiles and the TME to deepen our understanding of the potential applications of lncRNAs in liver cancer.

| Data sources
We retrieved TCGA genomic and transcriptomic data, including lncRNA data, for a total of 364 patients available in the TCGA-LIHC cohort. Multiple omics data and full clinical characteristics of the TCGA-LIHC cohort were downloaded from the Genomic Data Commons Data Portal. In addition, transcriptomic data and patient clinical information of the colon cancer (adenocarcinoma) cohort TCGA-COAD and stomach cancer (adenocarcinoma) cohort TCGA-STAD were obtained from TCGA and used for validation.

| Immune-related lncRNAs
Immune-related lncRNAs were obtained from the ImmLnc database (http://bio-bigda ta.hrbmu.edu.cn/ImmLnc). 9 The immunerelated lncRNAs were obtained in three steps. First, genome-wide and lncRNA-specific gene expression data were collected for the same patient. Then, all coding genes were ranked based on their correlation with each lncRNA expression. Finally, the genes in each immune-related pathway were mapped to a ranked gene list and the enrichment scores (ES) were calculated based on the Gene Set Enrichment Analysis (GSEA), which was converted to a lncRES score.
LncRNA-pathway pairs with a lncRES score of >0.995 and false discovery rate (FDR) of <0.05 were selected as immune-related lncR-NAs. The immune-related lncRNAs of liver cancer were downloaded from the ImmLnc database.

| Signature development of immunerelated lncRNAs
A total of 2908 immune-related lncRNAs in liver cancer were obtained, 699 of which were negatively correlated with immune-

| Molecular features
Single nucleotide variant and small insertion-deletion (INDEL) data were derived from the genomic data of the TCGA-LIHC cohort. The tumor mutation burden (TMB) and copy number variation (CNV) burden of the cohort were derived from data of a published study. 10 The top 10 genes with different somatic mutation types in the Lnc_ low and Lnc_high groups were displayed in a heatmap based on the frequency of the genes. The difference in gene mutation frequency between the two groups was compared using Fisher's exact test.
The genes were then grouped into eight known canonical pathways that included P53, WNT, RTK-RAS, TGFβ, Hippo, cell cycle, PI3K, and Notch, as previously described. 11 Samples in which genes of specific pathways contained somatic mutations were designated as having specific pathway alterations. Differences in the cancer pathway alteration frequency between the Lnc_low and Lnc_high groups were compared using Fisher's exact test.

| TME
Using genomic and transcriptomic data from the TCGA-LIHC cohort, we assessed the TME, including cancer immunograms, 28 immunecell subsets, five immune checkpoints, nine immune gene signatures, and immune subtypes. In accordance with a previous study, the cancer immunogram was designated using eight axes of the immunogram score (IGS), which reflect the seven steps in the cancer immunity cycle. 12 The IGS axes included IGS1, T-cell immunity; IGS2,

| Statistical analysis
The data were analyzed using R 3.6.1 and SPSS software version 24.0. Survival analysis was performed using Kaplan-Meier curves and compared between groups using a log-rank test. The chisquared test or Fisher's exact test was used to analyze the association between various genomic determinants. The student's t test was used to analyze the differences between the two groups when the data were normally distributed; otherwise, the Mann-Whitney U test was used. Statistical significance was set at p < 0.05.

| Immune-related lncRNA signature and prognosis in liver cancer
The significant lncRNA-pathway pairs in liver cancer with an absolute lncRES of >0.995 and an FDR of <0.05 were downloaded from the ImmLnc Database. A total of 2908 lncRNAs were significantly correlated with immune-related pathways, among which 699 were negatively correlated with immune-related pathways. We ranked these lncRNAs according to the sum of their ES in immune-related pathways. The 10 most significantly negatively correlated lncRNAs were selected and a lncRNA signature was obtained using logarith- the asscociation between etiology or liver fibrosis ishak score category and the lncRNA signature separately. However, we did not observe any significant difference among different subgroups ( Figure 2C,D).

| Distinct mutational profiles in the low and high lncRNA signature subgroups
To investigate the relationship between the lncRNA signature and mutational profiles, we analyzed the liver cancer genomic data of the TCGA-LIHC cohort. As shown in Figure 3A Figure 3B). Meanwhile, only five genes in the Lnc_low group had a higher number of mutations than that in the Lnc_high group. We also found that the number of mutations on genes involved in the p53 and WNT pathways was higher in the Lnc_high group compared with that in the Lnc_low group (two-sided Fisher's exact test, p < 0.01; Figure 3C). The other signaling pathways evaluated, including the RTK-RAS, TGFβ, Hippo, cell cycle, PI3K, and Notch pathways, were not significantly different between the two lncRNA signature groups (two-sided Fisher's exact test, p > 0.05; Figure 3C). Interestingly, we did not identify any significant difference in TMB and CNV burden between the lncRNA signature high and low subgroups ( Figure S2).

| Higher levels of the lncRNA signature correlated with an immune-excluded TME
We analyzed the transcriptomic data of the TCGA-LIHC data set to explore the association of the lncRNA signature with the immune microenvironment. First, to examine the immunological features between the high lncRNA signature group and low lncRNA signature group, we adopted a cancer immunogram that visually illustrates the state of the cancer immunity cycle. As shown in Figure 4A, four axes of the IGS, including IGS1, IGS3, IGS6, and IGS8, were significantly higher in the Lnc-low group compared with those in the Lnc-high group. The other four axes were not significantly different between the two groups. Next, we used single-sample GSEA (ssGSEA) to evaluate the relative abundance of 28 immune-cell subsets that infiltrated the tu-  Figure 4B). In addition, we observed a higher level of immunoinhibitory immune signature in the Lnc-high group ( Figure S3). We also analyzed the expression of immune checkpoint molecules between the two lncRNA groups and found that the expression levels of TNFRSF9, CD80, LAG3, PD-L1, and CTLA4 were significantly higher in the Lnc_high group ( Figure 4C). Finally, we obtained the immune subtype data of the TCGA-LIHC cohort analyzed by Thorsson et al. 10 and found a significant difference between the Lnc_high group and the Lnc_ low group (p < 0.01). The proportion of the C4 immune subtype was the highest in the Lnc_high group (51%), while the proportion of C3 immune subtype was the highest in the Lnc_low group (50%). The proportion of each immune subtype in the samples belonging to the two groups is shown in Figure 4D.

| DISCUSS ION
In recent years, significant progress has been made in the treatment of liver cancer, especially with the breakthrough of immunotherapy, which has allowed liver cancer intervention to enter a new era. Many advances have shown that the TME plays an important role in tumor development and is closely related to immunotherapy. 14,15 LncRNAs are emerging as critical regulators of the genome network and play critical roles in the development and activation of immune cells and immune regulation. 16 Here, we constructed a novel signature based on immune-related lncRNAs in patients with liver cancer that delineates an immune-excluded subtype of liver cancer with unfavorable clinical outcomes.
Immune-related lncRNA signatures have been explored in nonsmall cell lung cancer, breast cancer, and liver cancer. Meanwhile, co-expression analyses of lncRNAs and immune genes or immune cells have mainly been used to generate immune-related lncRNAs, showing that immune-related lncRNA signatures can be used to effectively predict the prognosis of patients. [17][18][19] Here, the immunerelated lncRNAs were obtained from the ImmLnc Database, which focuses on the lncRNAs that affect immune pathway activity. We constructed a lncRNA signature based on the top 10 inversely related immune-lncRNAs. The signature was shown to be correlated with poor prognosis (Figure 1), as well as progressed disease stage and histologic grade (Figure 2) in the TCGA-LIHC cohort data set.
The lncRNA expression was mainly correlated with the "cytokine receptors," "cytokines," "chemokines," and "antimicrobials" pathways, suggesting lncRNAs that regulate the immune pathway may play a crucial role in liver cancer progression and prognosis. Intriguingly, among the 10 lncRNAs, SLC16A1-AS1 have been reported correlated with poor prognosis in several cancers including liver cancer, 20,21 while others have not been reported yet. We validated the prognostic value of the lncRNA signature in colorectal and gastric cancer and found that the 10 lncRNA signatures were useful for predicting the prognosis of these cancers. This indicates that the lncRNA signature we generated from liver cancer data may have a universal function in gastrointestinal cancer.  23 Here, we found that the Lnc_high group had a significantly higher mutation frequency for both TP53 and CTNNB1 compared with that in the Lnc_low group. We also found that the number of mutations in the p53 and WNT pathways was higher in the Lnc-high group compared with that in the Lnc-low group.

Next-generation sequencing has provided a wealth of informa-
These findings suggest that high lncRNA signature levels may be associated with a cold immune microenvironment. TMB and CNV burden are well-known genomic signature related to immunological infiltrations, while we did not found any difference in TMB and CNV burden between the Lnc_high group and Lnc_low group. 24,25 This indicated that the associations between the lncRNA signature and immune microenvironment might be independent of TMB or CNV burden.
In recent years, with the development of checkpoint inhibitorbased immunotherapies, several immune checkpoint inhibitors have been applied in the clinical practice against liver cancer, including pembrolizumab, nivolumab, and atezolizumab, providing significant survival benefits to many patients. However, similar to non-small cell lung cancer and other cancers, only a small subset of patients F I G U R E 3 Distinct mutational profiles in low and high lncRNA signature subgroups. (A) Comparison of the top 10 frequently mutated genes between the Lnc_low group and Lnc_high group. (B) Comparison of the genes with significantly different mutation frequency between the Lnc_low and Lnc_high groups. (C) Comparison of eight canonical pathways between the Lnc_low and Lnc_high groups benefit from immune checkpoint blockade therapy. This is even true for combination immunotherapy, such as atezolizumab combined with bevacizumab, with less than one-third of patients in the Imbrave150 clinical trial showing a response. 26 Currently, commonly used biomarkers, including PD-L1 expression and tumor mutation burden, may play a role in patient selection for immunotherapy in some cancers, such as lung cancer, but they have very limited significance for liver cancer. The TME is increasingly recognized as a promising biomarker associated with immunotherapy responses. 5 Here, we found four axes of the IGS, IGS1, IGS3, IGS6, and IGS8, that were significantly higher in the Lnc_low group compared with those in the Lnc_high group. We also found that our immune-related lncRNA signature correlated with immune-cell infiltration. These findings suggest that immune-related lncRNA signatures are significantly correlated with the TME and may be integrated into a comprehensive biomarker system for immunotherapy.
In addition, we observed higher immunoinhibitor immune signature levels and expression levels of inhibitory checkpoints in the Lnc_high group compared with that in the Lnc_low group, which indicated that the tumors with Lnc_high signature were more likely to exhibit an immunosuppressive state and tended to escape immune cells within the TME. Recent studies have recognized that the overexpression of immune checkpoints, such as PD-L1, CTLA4, and LAG3, is associated with poor prognosis in liver cancer, 27,28 which is consistent with our result of the Lnc_high group having poor DFS and OS. Analyzing immune subtypes is another method to identify the composition of TME. The proportion of the C4 subtype was the highest immune subtype in the Lnc_high group, reaching 51%.
Meanwhile, the proportion of the C3 subtype was the highest immune subtype in the Lnc_low group, reaching 50%. This indicates that tumors of the Lnc_high group may represent an immuneexcluded subtype of liver cancer.
In conclusion, we constructed a novel immune-related lncRNA signature based on lncRNAs that affect immune pathway activity and determined that this lncRNA signature is correlated with the progression and prognosis of liver cancer. The immune-related ln-cRNA signature demonstrates a significant relationship with tumor mutations and the TME and may be integrated into a comprehensive biomarker system for immunotherapy. However, this immunerelated lncRNA signature needs to be applied to patients in a clinical setting to further verify its role in the prognosis of liver cancer and its ability to predict the efficacy of immunotherapy.

CO N FLI C T O F I NTE R E S T
The author(s) declared no potential conflict of interest with respect to the research, authorship, and/or publication of this article.

AUTH O R CO NTR I B UTI O N S
All listed authors participated in the study design, analysis, and interpretation of the data, drafted or revised the article, have agreed on the journal submitted currently, and given approval for the version to be published.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data can be obtained from TCGA database (https://portal.gdc. cancer.gov/).