Human endogenous retrovirus ERVK3-1 characterizes a metabolically active and immunosuppressive subtype of liver cancer.

Background
Liver hepatocellular carcinoma (LIHC) is one of the most common primary malignancies of the liver and a leading cause of cancer-related deaths worldwide, particularly in regions with high hepatitis B virus prevalence, such as China, Southeast Asia, and Africa [1]. Despite advancements in surgical resection, interventional treatments, targeted therapies (e.g., sorafenib, lenvatinib), and immune checkpoint inhibitors (e.g., anti-PD-1/PD-L1), the overall prognosis for LIHC patients remains poor, with a median overall survival of less than two years [2–4]. Increasing evidence suggests that tumor immune microenvironment (TIME) remodeling, metabolic dysregulation, and immune evasion mechanisms are central drivers of LIHC progression and therapeutic resistance [5].
Human endogenous retroviruses (HERVs) are sequences derived from ancient retroviral infections that have integrated into the human genome, accounting for approximately 8% of the human genome [6–8]. Among these, the HERV-K family is one of the most recent, structurally intact, and transcriptionally active subtypes [9, 10]. Aberrant activation of HERV-K has been implicated in the development of various human diseases and is emerging as a promising avenue for therapeutic intervention [11]. HERV activation can promote tumorigenesis through various mechanisms, including encoding tumor-associated antigens, inducing inflammatory responses, promoting oxidative stress, and interfering with cell cycle regulation [12]. Notably, HERVs also play a dual role in cancer immunity. While some HERV-encoded antigens may activate anti-tumor immune responses [13, 14], others—particularly envelope (ENV) proteins—can exert immunosuppressive effects that facilitate tumor immune escape. For instance, syncytin-2, an ENV protein, has been shown to suppress immune surveillance and promote tumor growth [15]. Additional immunosuppressive ENV proteins have been identified in other HERV families, including HERV-H, ERV3, HERV-P(b), and HERV-V [16].
ERVK3-1 (also known as HERVK_3p21.31), a member of the HERV-K family, has not yet been systematically studied in terms of its expression profile and functional role in tumors. Preliminary reports suggest that ERVK3-1 may be involved in transcriptional regulation and epigenetic changes in some cancers [17–19], but its expression pattern, functional role, and immunoregulatory functions in LIHC remain unexplored.
Therefore, this study systematically analyzes the expression levels of ERVK3-1 in LIHC and explores its associations with prognosis, metabolism, cell cycle, and the immune microenvironment. Through integrative analyses, we provide novel evidence that ERVK3-1 may contribute to LIHC progression by promoting immune evasion and T-cell exhaustion. These findings suggest that ERVK3-1 could serve as a potential biomarker and therapeutic target in LIHC.

Methods and materials

Data collection and processing
We retrieved the mRNA expression profiles and corresponding clinical information of LIHC patients from The Cancer Genome Atlas (TCGA) via the Genomic Data Commons (GDC) portal (https://portal.gdc.cancer.gov). Specifically, we downloaded the STAR-aligned read counts (HTSeq-Counts) and FPKM files from the GDC data release v33.0 (accessed on March 10, 2024). A total of 377 tumor samples were included, among which 50 cases had matched adjacent non-tumor tissues. Clinical features of these patients are detailed in Table 1. For further analysis, we retained 371 samples with both RNA-seq and clinical information. To normalize gene expression levels, we first converted raw counts to Transcripts Per Million (TPM) based on FPKM values and subsequently applied a log2(TPM + 1) transformation. This normalization was performed using R (version 4.2.3).

RNA-seq data processing
To assess the expression levels of ERVK3-1 and other HERV elements, we re-processed the RNA-seq data starting from raw FASTQ files using a customized pipeline designed to address the high sequence homology and frequent multi-mapping of HERV sequences. Quality control of raw reads was performed with FastQC (v0.11.9), followed by adapter trimming and removal of low-quality bases using Trimmomatic (v0.39). Cleaned reads were then aligned to the GRCh38 human reference genome with STAR aligner (v2.7.10a) in 2-pass mode to improve splice junction detection. To mitigate the issue of multi-mapped reads, we set the STAR parameter --outFilterMultimapNmax to 100 and retained multi-mapped reads with assigned mapping quality for subsequent quantification. Expression levels of both genes and repetitive elements, including HERVs, were quantified using TEtranscripts (v2.2.1), and the final expression matrices were normalized to TPM and log2(TPM + 1) for downstream analysis.

Single-cell RNA-seq data
In addition, we downloaded single-cell RNA sequencing datasets from the following sources: GSE16635 [20], GSE140228_10 × [21], and GSE146115 [22], all in .h5 format. The corresponding annotation results were obtained from TISCH (http://tisch.compbio.cn/home/). We used MAESTRO (v1.4.0) and Seurat (v4.3.0) for quality control, normalization, dimensionality reduction, clustering, and differential expression analysis. Cells were reclustered using the t-SNE algorithm, and ERVK3-1 expression was visualized across different annotated cell types.

Bioinformatic analyses

Differentially expressed genes (DEGs)
DEGs were identified using the Limma package (version 3.40.2) in R software. The results were visualized using the ggplot2 R package (version 3.3.6), and a volcano plot was generated to show the DEGs. A significance threshold of P < 0.05 and |Log2(Fold Change)| ≥ 1.5 (or greater than 1.0) was applied for DEG identification. Spearman’s correlation analysis was used to assess the relationship between ERVK3-1 expression and the top 10 DEGs.

Protein-protein interaction (PPI) network analysis and hub gene identification
Following DEG analysis, we constructed a protein-protein interaction (PPI) network using the online database STRING (https://www.string-db.org/). A confidence score threshold of 0.4 was set, while other parameters were kept at their default values. The resulting PPI network was visualized using Cytoscape software (version 3.9.1) [23]. Subsequently, we used the Cytoscape plugin Cytohubba with the MCC algorithm to identify the top 10 hub genes from the DEGs [24].

Functional enrichment analysis
To further investigate the potential functions of the target genes, we performed functional enrichment analysis. Gene Ontology (GO) is a widely used tool for annotating gene functions, specifically in terms of Molecular Function (MF), Biological Process (BP), and Cellular Component (CC). KEGG enrichment analysis is a practical method for analyzing gene functions and identifying high-level genomic functional information. To better understand the oncogenic role of the target genes, we utilized the ClusterProfiler package in R to analyze the GO functions of potential mRNAs and perform KEGG pathway enrichment analysis [25]. Additionally, we employed Gene Set Enrichment Analysis (GSEA) using the R package ClusterProfiler [version 4.2.1] [26]. We considered function or pathway terms with an adjusted P-value < 0.05 and false discovery rate (FDR) < 0.25 as statistically significant enrichment results.

Immune infiltration analysis
To evaluate the relationship between ERVK3-1 expression and immune cell infiltration in LIHC, we employed multiple complementary computational approaches: single-sample gene set enrichment analysis (ssGSEA), Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT), and The Tumor Immune Dysfunction and Exclusion (TIDE).
For ssGSEA, we quantified relative infiltration of 24 immune cell types across LIHC samples using the GSVA R package (v1.46.0) [27] and immune gene signatures from Bindea et al. [28]. For immune cell deconvolution, CIBERSORTx (https://cibersortx.stanford.edu/) was run with the LM22 signature matrix (22 immune cell types) [29], 100 permutations, and quantile normalization disabled (recommended for RNAseq). Immune cell fractions were computed using the CIBERSORT R script (v1.04). Tumor Immunophenotype (TIP) analysis (http://biocc.hrbmu.edu.cn/TIP/) integrated ssGSEA and CIBERSORT outputs to provide a comprehensive view of immune activity across the cancer immunity cycle.
We further examined correlations between ERVK3-1 expression and immunerelated genes (chemokines, interleukins, and immune checkpoint molecules) using Spearman’s correlation (stats package, base R v4.1.3). Immune response prediction was performed with the TIDE algorithm (http://tide.dfci.harvard.edu/) using default settings [30], and results were visualized with ggplot2 (v3.3.3) and ggpubr (v0.4.0).
Additionally, Immune Phenotype Scores (IPS) from The Cancer Immunome Atlas (TCIA) (https://tcia.at/) were used to estimate potential responsiveness to antiCTLA-4 and antiPD-1 therapies [31]. Enrichment analyses were carried out with clusterProfiler (v4.2.2), and only results with adjusted P < 0.05 were considered statistically significant.

Survival analysis
We used the Kaplan-Meier method to plot survival curves and assessed statistical significance using Cox regression analysis. The median expression level of ERVK3-1 was used as a cutoff value. Univariate and multivariate Cox regression analyses were performed to evaluate the impact of clinical variables on patient outcomes. In univariate Cox regression analysis, prognostic variables with a P-value < 0.05 were included in the multivariate Cox regression analysis. A forest plot was generated for visualization using the ggplot2 package.

Cell culture and experimental validation

Cell culture
THLE-2 cells were obtained from Fuheng Biotech Co., Ltd. (Shanghai, China) and cultured in Dulbecco’s Modified Eagle Medium (DMEM; Fuheng Biotech) supplemented with 10% fetal bovine serum (FBS; Tianhang Biotech Co., Ltd., Zhejiang, China) and 100 U/mL penicillin-streptomycin (Thermo Fisher Scientific, MA, USA). Huh-7 and HepG2 cells were purchased from the American Type Culture Collection (ATCC, Manassas, VA, USA) and maintained in DMEM (Invitrogen, Carlsbad, CA, USA) containing 10% FBS, 100 µg/mL penicillin, and 0.25 µg/mL streptomycin (Invitrogen). All cell lines were incubated at 37 °C in a humidified atmosphere containing 5% CO2.

siRNA transfection
Two siRNAs targeting ERVK3-1 and one scrambled negative control were synthesized by GenePharma (Shanghai, China) with the following sequences:

ERVK3-1 siRNA-1 (si-1): 5′-GGAGACCGGAGGGATCCAT-3′.

ERVK3-1 siRNA-2 (si-2): 5′-TCACTACTGGAGCGCCAGGCA-3′.

Scrambled siRNA (si-NC): 5′-UUCUCCGAACGUGUCACGUTT-3′.

Cells were transfected using Lipofectamine 2000 (Invitrogen) following the manufacturer’s instructions. Samples were collected 48 h after transfection for subsequent analyses.

RNA extraction and quantitative real-time PCR (qRT-PCR)
Total RNA was extracted using TRIzol reagent (Life Technologies, NY, USA) according to the manufacturer’s protocol. cDNA was synthesized with the PrimeScript RT Reagent Kit (Takara, Dalian, China). qRT-PCR was conducted using SYBR Premix Ex Taq (Takara) on a CFX96 Real-Time PCR System (Bio-Rad, CA, USA). Primer sequences were as follows: ERVK3-1: F 5’-TGGACAAACCGTGTGGGTGC-3’; R 5’-CCAGGAGTCCTGGGTTTTGC-3’ GAPDH: F 5’-AATGGGCAGCCGTTAGGAAA-3’; R 5’-GCGCCCAATACGACCAAATC-3’ PCR cycling conditions: 50 °C for 2 min, 95 °C for 2 min, followed by 40 cycles of 95 °C for 15 s and 60 °C for 1 min, with a final extension at 72 °C for 10 min. Relative gene expression was calculated using the 2^-ΔΔCt method, normalized to GAPDH.

Colony formation assay
After siRNA transfection, cells were seeded into 6-well plates at 500 cells/well (three technical replicates per group) and cultured for approximately 2 weeks with medium changed every 3 days. Colonies were washed with PBS, fixed with 4% paraformaldehyde for 20 min, and stained with 0.1% crystal violet. Images were acquired and colonies counted using a Cytation 5 imaging reader (BioTek Instruments, USA). Colony formation efficiency was calculated as the ratio of the number of colonies to the number of seeded cells.

Cell proliferation assay (CCK-8)
Transfected cells were seeded in 96-well plates at 1 × 10⁴ cells/well in six replicates, with PBS added to edge wells to prevent evaporation and blank wells containing only medium as controls. Starting from day 0, 10 µL of CCK-8 solution (Solarbio, Beijing, China) was added to each well daily, followed by incubation at 37 °C for 2 h. Absorbance was measured at 450 nm using a microplate reader, and corrected OD values were obtained by subtracting blank readings. Growth curves were plotted using OD values over 5 consecutive days to assess cell proliferation.

Statistical analysis
All analyses were performed using R software (v4.3.1). Differences in ERVK3-1 expression between tumor and adjacent normal tissues were evaluated by Wilcoxon rank-sum test (unpaired) or paired t-test, and correlations with clinical features, immune cell fractions, or gene expression levels were assessed using Spearman’s correlation. DEGs were identified with DESeq2, with P-values adjusted by Benjamini–Hochberg false discovery rate (FDR). Functional enrichment and immune-related analyses considered results with adjusted P < 0.05 as significant. For survival analysis, patients were stratified into high- and low-expression groups using the median ERVK3-1 value, an unbiased and widely used approach when no clinically validated cutoff is available. Kaplan–Meier curves with log-rank tests and Cox proportional hazards models were applied to estimate hazard ratios.
For experimental assays (qPCR, colony formation, and CCK-8), results are shown as mean ± standard deviation (SD) of at least three independent replicates, and statistical significance was determined by two-tailed Student’s t-test or one-way ANOVA (P < 0.05).

Results

Patient characteristics
We obtained clinical and RNA sequencing data for 377 LIHC patients from TCGA dataset. After excluding cases lacking clinical information or containing duplicate pathology records, 371 patients were included in the final analysis. Among them, 50 had matched adjacent normal tissue samples. Detailed clinicopathological characteristics are presented in Supplementary Table S1. To enhance the number of normal controls, we additionally incorporated gene expression data from 110 normal liver tissues in the GTEx dataset.

ERVK3-1 expression and prognostic relevance in LIHC
Pan-cancer analysis revealed that ERVK3-1 was significantly overexpressed in most cancer types. Specifically, expression in LIHC was markedly elevated compared to normal liver tissues, both in unpaired comparison (Fig. 1A) and in matched tumor–adjacent normal pairs (Fig. 1B). Focusing on LIHC, Cox proportional hazards regression demonstrated a significant association between high ERVK3-1 expression and reduced overall survival, suggesting its potential prognostic value (Fig. 1C). In parallel, receiver operating characteristic (ROC) curve analysis showed strong diagnostic performance in distinguishing LIHC from normal liver tissues, with an area under the curve (AUC) of 0.925 (95% CI: 0.889–0.961; Fig. 1D). To further validate these findings and address the limited number of normal liver samples in TCGA, we integrated gene expression data from the GTEx project. Consistent with the TCGA results, ERVK3-1 expression remained significantly elevated in LIHC tissues (P = 1.1e-43) (Fig. 1E), with a similar pattern observed in matched samples (P = 2.9e-13) (Fig. 1F). Consistent with the in silico analyses from TCGA and GTEx, in vitro validation further confirmed that ERVK3-1 mRNA expression was markedly elevated in hepatocellular carcinoma cell lines (Huh7 and HepG2) compared with the normal hepatic cell line THLE2 (Fig. 1G).

Association between ERVK3-1 expression and clinical characteristics in LIHC
We analyzed the correlation between ERVK3-1 expression levels and various clinical characteristics in LIHC samples from the TCGA database (Table 1). As shown in Fig. 2, ERVK3-1 was consistently overexpressed in various subgroups of LIHC patients compared to normal tissues, including by age, gender, race, BMI, T stage, N stage, M stage, pathological stage, histologic grade, AFP, albumin, prothrombin time, Child-Pugh grade, vascular invasion, residual tumor, and adjacent hepatic tissue inflammation.

To further delineate these associations, we performed univariate logistic regression analysis, which identified several clinical features significantly associated with ERVK3-1 expression. Specifically, lower expression was observed in patients > 60 years old (OR = 0.555, 95% CI: 0.367–0.838, P = 0.005), White individuals (vs. Asian and Black/African American; OR = 0.448, 95% CI: 0.294–0.683, P < 0.001), and those with BMI > 25 (OR = 0.502, 95% CI: 0.325–0.776, P = 0.002). Conversely, higher ERVK3-1 expression was significantly associated with more advanced disease features, including pathological T stage (T3&T4 vs. T1&T2; OR = 1.818, 95% CI: 1.125–2.938, P = 0.015), pathological stage (Stage III&IV vs. I&II; OR = 1.882, 95% CI: 1.149–3.083, P = 0.012), histologic grade (G3&G4 vs. G1&G2; OR = 3.163, 95% CI: 2.021–4.950, P < 0.001), and elevated AFP levels (> 400 vs. ≤400 ng/mL; OR = 2.756, 95% CI: 1.533–4.953, P < 0.001) (Table 2).

Prognostic value of ERVK3-1 in LIHC
To evaluate the prognostic relevance of ERVK3-1 in LIHC, we performed Kaplan–Meier survival analysis based on its expression levels (Figs. 3A-C). Patients were stratified into high and low expression groups using the median ERVK3-1 expression as the cutoff. The analysis revealed that high ERVK3-1 expression was significantly associated with poorer overall survival (OS) (hazard ratio [HR] = 1.52, 95% CI: 1.07–2.15, P = 0.018) and disease-specific survival (DSS) (HR = 1.56, 95% CI: 1.00–2.44, P = 0.047). No significant association was observed with progression-free interval (PFI).

To further investigate independent prognostic factors, we conducted univariate and multivariate Cox regression analyses. In the univariate model, pathologic T3 stage (adjusted HR = 2.674, 95% CI = 1.761–4.060, P < 0.001), T4 stage (adjusted HR = 5.386, 95% CI = 2.690–10.784, P < 0.001), M stage (adjusted HR = 4.077, 95% CI = 1.281–12.973, P = 0.017), residual tumor (R2, adjusted HR = 11.749, 95% CI = 1.595–86.516, P = 0.016), and ERVK3-1 (adjusted HR = 1.524, 95% CI = 1.076–2.159, P = 0.018) were significantly associated with worse OS. Multivariate analysis confirmed that pathologic T3 stage (adjusted HR = 3.044, 95% CI = 1.824–5.079, P < 0.001), T4 stage (adjusted HR = 5.168, 95% CI = 1.884–14.174, P = 0.001), and high ERVK3-1 expression (adjusted HR = 1.739, 95% CI = 1.114–2.715, P = 0.015) were independent predictors of poor OS in LIHC (Fig. 3D).

Construction and validation of a prognostic nomogram
To facilitate individualized survival prediction, we constructed a prognostic nomogram integrating the independent risk factors—pathologic T stage and ERVK3-1 expression—identified in multivariate analysis. The nomogram assigns a weighted score to each variable, with the total score corresponding to a predicted overall survival probability, where higher scores indicate poorer prognosis (Fig. 4A).The predictive performance of the nomogram was evaluated using the bootstrap-corrected concordance index (C-index), which yielded a value of 0.644 (95% CI: 0.617–0.672), indicating moderate accuracy. Calibration plots demonstrated strong concordance between predicted and observed survival at 1-, 3-, and 5-year intervals (Figs. 4B–D), supporting the model’s reliability. Furthermore, time-dependent ROC curve analysis showed that the combination of T stage and ERVK3-1 expression achieved reasonable discriminative power in predicting survival outcomes (Fig. 4E). Collectively, these results suggest that the proposed nomogram provides a robust and practical tool for prognostic assessment in LIHC.

Identification of differentially expressed genes (DEGs) in LIHC and PPI network analysis
To elucidate the molecular pathways associated with ERVK3-1 in LIHC, we stratified patients into high and low expression groups based on the median ERVK3-1 level and conducted differential gene expression analysis. A total of 1,316 DEGs were identified, including 1,010 downregulated genes (76.7%) and 306 upregulated genes (23.3%) (adjusted p-value < 0.05, |Log2-FC| >1.5) (Fig. 5A). A heatmap of the top 50 upregulated and 50 downregulated genes is shown in Fig. 5B. The top 10 DEGs most strongly correlated with ERVK3-1 expression were TSPAN6, TNMD, DPM1, SCYL3, C1orf112, FGR, CFH, FUCA2, GCLC, and NFYA (Fig. 5C, Supplementary Table S2). To further explore the functional relevance of these genes, we constructed a protein-protein interaction (PPI) network using the STRING database (Supplementary Table S3). Using the MCC algorithm in Cytoscape, the top 10 hub genes identified were: CCNB1, CDC20, CCNA2, KIF11, CDK1, BUB1, DLGAP5, KIF2C, BUB1B, and KIF20A (Fig. 5D). These hub genes may represent key nodes in ERVK3-1–related regulatory networks in LIHC.

Functional enrichment analysis: GO, KEGG, and GSEA
To explore the potential biological roles of ERVK3‑1 in LIHC, we performed functional enrichment analyses based on differentially expressed genes. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses showed that genes upregulated with high ERVK3‑1 expression were enriched in small molecule catabolism, oxidoreductase activity, and mitochondrial matrix, whereas downregulated genes were mainly involved in cell cycle regulation and DNA replication (Fig. 6A–D; Supplementary Table S4). Consistent with these findings, correlation analysis revealed that ERVK3‑1 expression was positively associated with key regulators of cell cycle checkpoints and the DNA damage response, including ATM, ATR, CHEK1/2, TP53, and multiple DNA repair effectors (P < 0.05; Supplementary Fig. 1). This association suggests that ERVK3‑1 may participate in metabolic reprogramming and cell‑cycle control under DNA damage stress, potentially conferring adaptive advantages to tumor cells facing genomic instability. Supporting these observations, Gene Set Enrichment Analysis (GSEA) further highlighted pathways related to amino acid catabolism, oxidative metabolism, and P450‑mediated drug metabolism (Fig. 6E–F; Supplementary Table S5).

To experimentally validate these implications, ERVK3‑1 was silenced in Huh7 and HepG2 cells using siRNA, followed by colony formation and CCK‑8 assays. Knockdown markedly reduced colony formation in Huh7 cells (24.33% vs. 57.33%, P < 0.05) and significantly slowed proliferation (Supplementary Fig. 2A–D). A similar trend was observed in HepG2 cells, with colony formation decreasing from 87.66 to 63.44% (P < 0.05) and growth inhibition confirmed by CCK‑8 analysis (Supplementary Fig. 2E–F).

Correlation between ERVK3-1 expression and immune infiltration in LIHC
Given the emerging role of human endogenous retroviruses (HERVs) in tumor immunity, we investigated the association between ERVK3-1 expression and immune cell infiltration in LIHC. Single-sample Gene Set Enrichment Analysis (ssGSEA) revealed that ERVK3-1 expression was negatively correlated with several immune cell populations, including Th17 cells (r = -0.289, P < 0.001), dendritic cells (DC, r = -0.236, P < 0.001), neutrophils (r = -0.230, P < 0.001), γδ T cells (Tgd, r = -0.154, P < 0.01), regulatory T cells (Tregs, r = -0.136, P < 0.01), mast cells (r = -0.117, P < 0.05), and NK cells (r = -0.108, P < 0.05). In contrast, positive correlations were observed with Th2 cells (r = 0.293, P < 0.001), NK CD56bright cells (r = 0.213, P < 0.001), follicular helper T cells (TFH, r = 0.196, P < 0.001), T helper cells (r = 0.156, P < 0.01), and activated dendritic cells (aDCs, r = 0.113, P < 0.05) (Fig. 7A).

These findings were corroborated by CIBERSORT analysis, which identified negative correlations between ERVK3-1 expression and mast cells resting (r = -0.225, P < 0.001), monocytes (r = -0.137, P < 0.01), macrophages M1 (r = -0.129, P < 0.05), and M2 (r = -0.117, P < 0.05). Conversely, positive correlations were found with follicular helper T cells (r = 0.225, P < 0.001), regulatory T cells (r = 0.161, P < 0.01), and CD4 memory activated T cells (r = 0.144, P < 0.01) (Fig. 7B). To further elucidate ERVK3-1’s role in modulating the immune microenvironment, we analyzed its correlation with chemokines and cytokines. ERVK3-1 expression showed significant associations with multiple chemokine and interleukin cytokine genes, suggesting its broad involvement in immune signaling pathways within LIHC (Figs. 7C–D). Finally, we examined the relationship between ERVK3-1 and the cancer-immunity cycle using the Tracking Tumor Immunophenotype (TIP) database. ERVK3-1 expression was positively associated with early steps (steps 1–4) of the cycle, including antigen release, presentation, and immune cell recruitment. However, it was negatively correlated with later steps, such as immune cell infiltration into tumors, T-cell recognition of cancer cells, and tumor cell killing (Fig. 7E).
Collectively, these results indicate that ERVK3-1 is intricately linked to immune cell dynamics and may contribute to immune evasion and poor clinical outcomes in LIHC.

Potential role of ERVK3-1 in immune evasion and immunotherapy response in LIHC
To explore the potential mechanisms underlying ERVK3-1’s role in immune evasion, we analyzed single-cell RNA-seq data from public databases and identified high expression of ERVK3-1 in immune subpopulations such as T proliferative cells (Tprolif), fibroblasts, and CD8 + exhausted T cells (CD8Tex). Among these, ERVK3-1 expression was most pronounced in the Tprolif subpopulation across multiple datasets, including GSE166635, GSE140228_10x, and GSE146115 (Fig. 8).

Moreover, ERVK3-1 expression was positively correlated with several inhibitory immune checkpoint molecules, including CD160, CD96, CSF1R, CTLA4, HAVCA2, IL10RB, LAG3, LGALS9, TIGIT, and VTCN1 (Fig. 9A), suggesting a possible role for ERVK3-1 in mediating T-cell suppression and immune evasion. Consistently, ERVK3-1 expression also showed positive associations with key markers of T-cell exhaustion, including PDCD1 and HAVCR2, as well as with effector molecules such as GZMB, GZMK, PRF1, and IFNG (Fig. 9B). This expression profile suggests that ERVK3-1 may contribute to the establishment of an exhausted T-cell phenotype, thereby weakening anti-tumor immune responses and facilitating immune escape.

To evaluate the clinical relevance of ERVK3-1 in predicting immunotherapy response, we calculated Tumor Immune Dysfunction and Exclusion (TIDE) scores. The high ERVK3-1 expression group exhibited significantly elevated TIDE scores (P = 1e-03), indicating a higher likelihood of immune evasion and a lower probability of benefiting from immune checkpoint inhibitor (ICI) therapy (Fig. 9C). Additionally, based on data from The Cancer Immunome Atlas (TCIA), we analyzed the Immune Phenotype Score (IPS) to predict responsiveness to anti-CTLA-4 and anti-PD-1 therapies. In the non-responder group, patients with low ERVK3-1 expression had a lower likelihood of resistance to ICIs (P = 0.01) (Fig. 9D). However, among the responder group, no significant difference in response was observed based on ERVK3-1 expression levels.
Collectively, these findings suggest that high ERVK3-1 expression may contribute to immune evasion by promoting T-cell exhaustion and upregulating immune checkpoint pathways. Moreover, ERVK3-1 may serve as a potential predictive biomarker for identifying LIHC patients who are less likely to benefit from ICI-based immunotherapy.

Discussion
LIHC ranks as the third leading cause of cancer-related mortality worldwide, characterized by high recurrence and metastatic potential that severely impairs patient prognosis [32]. Recent advances in immune checkpoint blockade (ICB) therapies, such as nivolumab (anti-PD-1) and atezolizumab (anti-PD-L1), have transformed treatment strategies for advanced LIHC [33]. Despite this, primary resistance to these therapies is common, underscoring the need for a deeper understanding of the tumor immune microenvironment, particularly T-cell exhaustion.
As relics of ancestral viral insertions, human endogenous retroviruses (HERVs) are now recognized for their capacity to modulate tumor immunogenicity and shape immune responses, highlighting their potential significance in LIHC pathogenesis [34, 35]. Although several HERV families—including HERV-K, HERV-H, and ERV3—have been implicated in tumor progression via mechanisms such as immune evasion, chronic inflammation, and epigenetic alterations [12, 36], their specific roles in LIHC remain underexplored. In this study, we systematically characterized ERVK3-1, an intronic element of the HERV-K family, and found it to be significantly associated with metabolic reprogramming, cell cycle activation, and T cell exhaustion in LIHC. These functional patterns differ from traditional HERV protein-mediated effects, suggesting an alternative, non-protein-mediated role of ERVK3-1 through host gene regulation [37]. Importantly, ERVK3-1 displayed robust diagnostic performance (AUC = 0.925) and demonstrated independent prognostic value, underscoring its potential as both a biomarker and therapeutic target in LIHC [38]. Consistent with findings of HERV-K overexpression across multiple malignancies [39–46], our results extend previous LIHC-specific studies such as Ma et al. [38] by identifying ERVK3-1 as a distinct and prognostic biomarker.
HERV-K elements have been implicated in various oncogenic processes. In pancreatic cancer, silencing of HERV-K env RNA suppresses tumor growth via inhibition of the Ras-ERK-RSK pathway [43]. Similarly, HERV-K-targeted CAR-T cells reduce breast cancer proliferation and metastasis by altering p53, MDM2, and ERK signaling [47]. In line with these findings, our DEGs analysis identified key regulators of the cell cycle and mitosis, such as CCNB1, CCNA2, CDK1, KIF11, KIF2C, and KIF20A, all of which are implicated in chromosomal segregation, spindle formation, and cell cycle progression [48–51]. High ERVK3-1 expression may activate these pathways, promoting rapid tumor cell division and proliferation. Moreover, our analyses in GO and KEGG further revealed that low ERVK3-1 expression is associated with cell cycle regulation, while high expression is linked to metabolic pathways, including oxidative-reductase activity, mitochondrial function, and cytochrome P450 metabolism. Mitochondrial dysfunction and oxidative stress are known contributors to LIHC progression, with the hypoxic tumor microenvironment triggering mitochondrial fragmentation and immune escape mechanisms [52, 53]. P450 enzymes also influence tumor drug resistance and immune evasion through metabolic remodeling [54, 55].
These observations point to a potential coupling between metabolic reprogramming and immune evasion. Metabolic alterations, such as redox imbalance and P450 activity, may dysregulate CDK1 and cell cycle proteins, contributing to immune evasion. This metabolic-cell cycle imbalance facilitates immune suppression via PD-L1 upregulation and ROS/lactate accumulation [56]. Therefore, we hypothesize that ERVK3-1 may influence tumor cell growth through the regulation of mitochondrial function, metabolic reprogramming, and the cell cycle protein network. It could also impact the expression of immune checkpoint molecules and the remodeling of the tumor immune microenvironment (TME) via a metabolic-cell cycle-immune coupling mechanism. This suggests that ERVK3-1 may influence tumor progression by modulating mitochondrial function and metabolic reprogramming, contributing to immune evasion.
Notably, high ERVK3-1 expression correlates with increased levels of Tregs and various chemokines (e.g., CCL20, CCL22, CXCL16) and immune-regulatory cytokines (e.g., IL-10, IL-18), suggesting its role in immune suppression in liver cancer [57, 58]. HERVs can promote immune tolerance via TERT-mediated IFN signaling and Treg recruitment [59–61]. In retrovirus-infected mouse models, FoxP3⁺ Tregs upregulate IL-10, contributing to persistent immunosuppression [62]. Subsequent analysis using the TIP tool to assess its function across the seven stages of the cancer immunity cycle revealed a “pro-immune then suppressive” dynamic [63]. This shift in immune response suggests that ERVK3-1 may regulate the temporal and spatial organization of the immune microenvironment, thereby promoting the establishment of immune tolerance [64]. These findings linking ERVK3-1 expression with the upregulation of immunosuppressive factors such as Tregs and IL-10 further support our observation that ERVK3-1 is involved in regulating immune tolerance and poor prognosis in LIHC patients, positioning it as a potential therapeutic target for immune intervention.
Single-cell RNA-seq data further confirmed ERVK3-1 enrichment in immunosuppressive subpopulations, including Tprolif and CD8⁺ exhausted T cells. Its expression correlates with exhaustion markers and effector genes, indicating persistent antigen-driven T-cell dysfunction [65]. This quantitative and spatial shift from “activation to exhaustion” could represent a key mechanism underlying ERVK3-1’s role in immune escape. TIDE and IPS analyses also indicated that ERVK3-1 overexpression is associated with lower immunotherapy responsiveness and worse prognosis.
Taken together, our findings suggest that ERVK3-1 plays a multifaceted role in LIHC progression through the regulation of metabolism, cell cycle dynamics, and immune evasion. Despite these insights, the study is limited by its reliance on computational analyses without experimental validation. Nevertheless, our work integrates bulk RNA-seq, single-cell transcriptomics, and targeted in vitro experiments, providing a coherent and internally consistent framework that robustly supports our conclusions. Future studies will further employ genome editing technologies such as CRISPR/Cas9 or RNA interference to directly interrogate the causal role of ERVK3‑1 and clarify whether it functions as a driver or a passenger in hepatocarcinogenesis. Complementary assays, including metabolic characterization and immune co-culture models, will also be used to define its mechanistic impact and therapeutic potential.
ERVK3-1 demonstrates potential as a diagnostic and prognostic biomarker and may guide immunotherapy strategies in LIHC. While promising, clinical translation of ERVK3-1 requires standardized detection and broader validation. Further mechanistic insights are also necessary to develop effective targeted interventions. Clinical investigations are warranted to establish ERVK3-1 as a biomarker and therapeutic target in immunotherapy settings. Overall, ERVK3-1 emerges as a promising candidate for precision oncology in liver cancer, bridging metabolic regulation and immune modulation within the tumor microenvironment.

Conclusion
In summary, this study presents a comprehensive transcriptomic and immunogenomic analysis of ERVK3-1 in LIHC. Our findings suggest that ERVK3-1 may promote tumor progression by regulating key biological processes such as metabolism, the cell cycle, and immune evasion. Its strong association with immunosuppressive markers and T-cell exhaustion highlights its role in shaping an immunosuppressive tumor microenvironment. Furthermore, high ERVK3-1 expression correlates with higher TIDE scores and lower immunotherapy responsiveness, indicating its potential as a biomarker for poor prognosis and resistance to immune checkpoint blockade. As a member of the HERV-K family, ERVK3-1 may represent a novel link between metabolic and immune dysregulation in LIHC, offering new avenues for therapeutic targeting. Further experimental validation and mechanistic studies are warranted to clarify its clinical utility.

Supplementary Information
Below is the link to the electronic supplementary material.