The mRNA transcriptional data of HCC and normal liver tissues were obtained from The Cancer Genome Atlas (TCGA) (http://gdc.cancer.gov) datasets. All raw data were further analyzed after being standardized by log₂ transformation. The R (version 4.0.3) and GraphPad Prism (version 7.0) software were utilized for analysis. We chose |log₂ fold change (FC)|≥1 and adjusted p-value <0.05 as statistically significant genes.

A total of 113 HCC patients who underwent curative surgery in Eastern Hepatobiliary Surgery Hospital affiliated to Second Military Medical University from 2008 to 2014 were enrolled in this study. They were followed up postoperatively until December 2020. Patients confirmed of HCC by histopathological diagnosis and without adjuvant anticancer treatment such as radiotherapy and chemotherapy before surgery were chosen in the study. The clinicopathological data of selected patients were gained from the patient’s hospitalization records. The studies have obtained ethics approval by medical ethics committees of Tongji Hospital and Eastern Hepatobiliary Surgery Hospital, and all the patients involved in the study signed an informed consent form.

Total RNA in liver of HCC patients and paired adjacent normal tissues was extracted by TRIzol reagent (Sigma-Aldrich, USA) under the instruction of the manufacturer’s protocol. The cDNAs were synthesized using a Prime Script TM RT kit (Takara, Japan) in reverse transcription reactions and the results were normalized to endogenous glyceraldehyde-3-phosphate dehydrogenase (GAPDH) expression. The primer sequences of ELOVL1 are the following: FORWARD (5′-3′) TTATTCTCCGAAAGAAAGACGGG, REVERSE (5′-3′) TTATTCTCCGAAAGAAAGACGGG. The primer sequences of ELOVL3 are the following: FORWARD (5′-3′) CTGTTCCAGCCCTATAACTTCG, REVERSE (5′-3′) GAATGAGGTTGCCCAATACTCC. The primer sequences of GAPDH are the following: FORWARD (5′-3′) AATGGGCAGCCGTTAGGAAA, REVERSE (5′-3′) GCGCCCAATACGA CCAAATC.

Total protein was extracted with RIPA lysis buffer (Beyotime, China) mixed with a protease inhibitor (Beyotime, China) and centrifuged at 14,000g for 15 min at 4°C. Supernatant were quantified by the BCA Protein Assay Kit (Beyotime, China). Equal amounts of protein were separated, transferred onto polyvinylidene difluoride membranes, then blocked with 5% bovine serum albumin for 2 h at room temperature and incubated with primary antibodies against ELOVL1 (diluted 1:1,000, TA0670, Abmart, China) and GAPDH (diluted 1:25,000, 60004-1, Proteintech). The membrane was further reacted with horseradish peroxidase–conjugated secondary antibody (diluted 1:5,000, Beyotime, China) for 1 h at room temperature. The band intensity was analyzed by the ChemiDoc XRS systems (Bio-Rad Laboratories, United States) and Image J software.

The tissue microarray was constructed using the paraffin-embedded tissue of the 113 patients. Tumor sections were incubated with primary anti-ELOVL1 (diluted 1:200, TA0670, Abmart, China) antibody in fridge overnight at 4°C. Then, the tissues were incubated with the secondary antibody (1:1,000 dilution, Thermo Fisher Scientific, A-10042, Massachusetts, USA) at 37°C for 1 h, and then covered by 3,3-diaminobenzidine (ZLI-9032, Zhongshan Biotech, China). Subsequently, all tissues were reviewed using the light microscope (Olympus 600 auto-biochemical analyzer, Japan). The positive cells score negative: 0%–5%; low: 6%–25%; medium: 26%–50%; high: >50%.

To determine the differentially expressed gene (DEG) pattern between HCC patients with high and low ELOVL1expression, HCC patients in TCGA database were divided into high and low groups in line with the median ELOVL1 expression value. DEGs were determined by using the Limma package with the absolute value of logFC (logFoldchange) ≥1 and p-value ≤0.05. The functions of these identified DEGs were explored by hallmark gene sets, i.e., Gene Ontology (GO) gene sets and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Gene set enrichment analysis (GSEA) was also analyzed using “clusterProfiler” package in R. The results were visualized using cluster Profiler and ggplot2 R packages (threshold: P < 0.05).

CIBERSOR was applied to explore the expression of ELOVL1 and relationship with the abundance of 22 tumor-infiltrating immune cells (TILCs) including CD8⁺ T cells, CD4⁺ T cells, B cells, neutrophils, macrophages, and dendritic cell by gene expression profiling in 44 cancer types from 10,180 samples (18). In addition, the TIMER database (https://cistrome.shinyapps.io/timer/) was also applied to inspect the correction between ELOVL1 and the abundance of TILCs through the “gene” module (19). The mRNA level of ELOVL1 and relationship with 60 immune checkpoints (24 inhibitory and 36 stimulatory) was identified by TISIDB database (http://cis.hku.hk/TISIDB/index.php) (20). The tumor immune dysfunction and exclusion (TIDE) was conducted to explore the correlation of ELOVL1 with TIDE scores and to predict the response possibility between ELOVL1-low and ELOVL1-high group to immune checkpoint inhibitor (21). Additionally, we analyzed the relationships of ELOVL1 expression with several clinical cohorts received immunotherapy through the TIDE website.

To identify the regulations of ELOVL1 expression and HCC, we first integrated ELOVL1 gene expression in TCGA pan-cancer (PanCAN, n = 10,535, g = 60,499). Second, as described in previous articles, TMB, MSI, homologous recombination deficiency (HRD), and neoantigen were used to evaluate the relationship between ELOVL1 and tumor mutation and treatment sensitivity (22). For TMB, we focused on the relationship between high expression of ELOVL1 and TMB. Last, mutant-allele tumor heterogeneity (MATH), purity, ploidy, and loss of heterozygosity (LOH) were used to inspect the association between ELOVL1 and tumor heterogeneity, as described in previous articles (23).

Survival curves were generated from the data in TCGA. The correlations between ELOVL1 expression and clinicopathological features were evaluated by Chi-square test evaluation. Kaplan–Meier survival curves were constructed using R (version 4.0.3). P < 0.05 was thought about statistically significant in the current study.

At pan-cancer level, TCGA and GTEx database were utilized to explore the expression of ELOVLs between tumor and normal tissues. The results showed that the mRNA level of ELOVLs in most tumor tissues including HCC was higher than that in normal tissues ( Figures 1A–G ). Compared with normal tissues, ELOVL1, ELOVL2, ELOVL3, ELOVL5, and ELOVL7 were significantly increased in HCC (left subgroups, Figures 2A–G ). In addition, ELOVL1, ELOVL2, ELOVL3, ELOVL4, and ELOVL5 were significantly overexpressed in the paired comparison between HCC and adjacent tissue (right subgroups, Figures 2A–G , Figure 2H ). These data suggested that the ELOVLs family might play an important role in HCC.

Furthermore, univariate Cox regression ( Figure 3A ) and multivariate Cox regression ( Figure 3B ) analyses were used to analyze the role of ELOVLs family on the survival outcomes of HCC. Among the ELOVLs family, only ELVOL1 and ELOVL3 might have an impact on the survival of liver cancer patients. Consistent with this, the high expression of both ELVOL1 and ELOVL3 would reduce the overall survival (OS) ( Figures 3C, E ) and disease-free survival (DFS) ( Figures 3D, F ) of patients with HCC. These results suggested that ELVOL1 and ELOVL3 in the ELOVLs family were strongly associated with HCC and have the potential to predict the prognosis of patients with HCC.

To validate the results of the bioinformatics analysis in the ELOVLs family, we collected human samples of liver cancer and adjacent tissues. The results of RT-qPCR showed that the mRNA expression of ELOVL1, but not ELOVL3, was higher in HCC patients than of paired adjacent normal liver tissues ( Figures 4A, B ). Consistent with this, Western blot results showed that the protein expression of ELOVL1 was higher in HCC patients than of paired adjacent normal liver tissues ( Figures 4C, D ).

Furthermore, the immunohistochemical chip of 113 liver cancer samples showed 68 cases of them displayed high expression of ELOVL1 levels and 45 cases displayed low ELOVL1 expression levels ( Figures 5A ). Combined with bioinformatics analysis and human samples validation, these data indicated that ELOVL1 was significantly upregulated in HCC tissues compared with adjacent normal liver tissues, which suggested that ELOVL1 might play an important role in HCC.

Based on the expression profile of ELOVL1, we used transcriptome data from TCGA database to reclassify HCC patients into high-expressed ELOVL1 and low-expressed ELOVL1 groups. Then, logistic regression model was constructed to explore the relationship of ELOVL1 and clinicopathological characteristics in HCC patients. The results showed that ELOVL1was strongly associated with the tumor grade and tumor T stage of HCC patients ( Table 1 ). Then, Cox proportional hazard models were constructed to explore the effects of ELOVL1 on patient survival times. The univariate Cox analysis showed that tumor stage (HR = 1.144, p < 0.001), T stage (HR = 1.219, p < 0.001), metastasis (HR = 1.273, p =0.010), N (regional lymph node) (HR = 1.232, p =0.024), and ELOVL1 expression (HR = 1.266, p < 0.001) were independent factors for OS of HCC patients ( Figure 5B ). Multivariate Cox regression analysis showed that ELOVL1 expression was confirmed to be a statistically significant predictor of OS in HCC patients (HR = 1.011, p < 0.001; Figure 5C ). Highly expressed ELOVL1 in HCC patients would increase the survival risk and reduce the survival probability. In addition, T stage (HR = 1.199, p < 0.001) and metastasis (HR = 1.309, p =0.032) were also confirmed as an independent risk factor for OS.

In conclusion, in terms of clinical applications, ELOVL1 was closely related to tumor grade and tumor T stage and might predict the prognosis of HCC patients.

To identify how ELOVL1 participated in the tumorigenesis of HCC, we applied hallmark, GO, and GSEA analysis to explore the signaling pathways that involved in HCC. The hallmark enrichment analysis implied that ELOVL1 was involved in the E2F targets, G2M_checkpiont, mototic spindle, myc targets, MTORC1 signaling, DNA repair, unfolded DNA repair, and protein section. The GSEA study showed that ELOVL1 participated antigen processing and presentation (ES = 0.8987, NP = 0.0048), pathways in cancer (ES = 0.6627, NP=0.0124), leukocyte transendothelial migration (ES = 0.7697, NP = 0.0136), NOD-like receptor signaling pathway (ES = 0.8072, NP = 0.0187), and Toll-like receptor signaling pathway (ES = 0.7631, NP = 0.0260) ( Supplementary Figure 1 ). Most of these pathways were related to inflammation and immune response, which is closely correlated with the tumorigenesis of HCC. Therefore, we further evaluated the relationship between ELOVL1 and the immune-related pathways. The results showed that ELOVL1 participated the antigen processing and presentation and innate immune response ( Figure 6A ). Besides, ELOVL1 also participated in major histocompatibility complex protein binding and adaptive immune response ( Figures 6B–D ). The results implied that ELOVL1 be involved in the tumorigenesis of HCC through immune microenvironment.

Based on the results of functional enrichment, we analyzed the immune components between high expression ELOVL1 and low expression ELOVL1 groups, which included immune cell infiltration and immune checkpoints. The results showed the mRNA expression of ELOVL1 was closely related to the immune cell infiltration, including CD8⁺ T cells, CD4⁺ T cells, neutrophils, macrophages, and dendritic cells through both the CIBERSOR and TIMER ( Figures 7A–C and Supplementary Figure 2 ). In terms of immune checkpoints, the heat maps illustrated that ELOVL1 were closely connected with immune checkpoint. ELOVL1 was positively correlated with PD-1 (PDCD1), CTLA4, LAG3, endothelial growth factors (VEGFs), and so on, which are the common immune checkpoints in the HCC tumor-immune microenvironment (TIME) ( Figure 7D ).

Next, we performed the TIDE analysis of the effect of ELOVL1 on immunotherapy. These results showed that ELOVL1 had an obvious positive correlation with the TIDE score ( Figures 8A, B ). Specifically, patients with elevated ELOVL1 levels were more likely to show no response to immunotherapy according to the results from the TIDE analysis, whereas patients with low ELOVL1 level may be more sensitive to immunotherapy treatment ( Figure 8C–D ). All these results implied that ELOVL1 might participate in HCC via immune evasion and might predict the efficacy of immunotherapy.

.

Studies have shown that indicators of genome heterogeneity might be the biomarker of immunotherapy, so we analyzed the correlations between ELOVL 1 and genome heterogeneity such as MSI, TMB, and neoantigen in HCC (24–26). As shown in Figures 8E–H , ELOVL1 is positively related to MSI and HRD (R = 0.298967888875077, P = 1.71923786304186e-8) in eight cancers including HCC (R = 0.104455471439159, P = 0.0455298564418963). Meanwhile, ELOVL1 is negatively related to TMB (R = −0.12, P = 0.047). However, ELOVL1 was not connected with neoantigen. MSI and TMB are the most studied tumor immunotherapy markers in cancers including HCC. In the correlation between ELOVL1 and MSI, TMB illustrated that ELOVL1 might be the biomarker of the efficacy of immunotherapy, which is consistent with the former results. As for the tumor heterogeneity, ELOVL1 was positively related to MATH (R = 0.15, P = 0.020), tumor ploidy (R = 0.12, P = 0.029), and LOH (R = 0.38, P = 0.001) in HCC. Additionally, ELOVL1 was negatively related to tumor purity (R = −0.11, P = 0.035) in HCC ( Supplementary Figure 3 ). A few studies implied that oncogenic pathways driven by genetic mutations were also related to immune microenvironment, which, in turn, could impact response to immunotherapies in several types of cancer including HCC (27, 28). Thus, the mutation database was used to analysis the associations between ELOVL1 and mutational landscape in HCC. As shown in Supplementary Figure 4 , ELOVL1 was related to ADRA1D, POTEH, GRIN1, NLRP12, SHANK1, MYH7, ZFPM2, FAM47A, TSC2, RB1, CSMD1, CTN, NB1, and TP53, which were common mutations in HCC. The relationship of MSI, TMB, and mutational landscape indicated that ELVOL1 had the potential to predict the efficacy of immunotherapy of HCC.