Introduction

Clear-cell renal cell carcinoma (ccRCC), the most prevalent form of kidney cancer, accounts for a significant proportion of cancer-related fatalities.1 The well-documented somatic mutations in the VHL gene leading to activation of angiogenesis, along with recent discoveries highlighting dysregulation of tumor microenvironment (TME) have propelled the Food and Drug Administration (FDA)’s approval of various therapies for metastatic RCC.2 Recent phase III trials have informed current guidelines, which advocate for the initial treatment with either dual immune checkpoint inhibition (ICI) or a combination therapy comprising a vascular endothelial growth factor receptor tyrosine kinase inhibitor (TKI) and an anti-PD-1 ICI.3

The pivotal CheckMate-214 trial demonstrated a 42% response rate to the combination of ipilimumab and nivolumab (Ipi/Nivo). However, an early disease progression was observed in nearly 20% of patients.4 Despite efforts, including advanced single-cell RNA sequencing (RNA-seq), to identify predictive biomarkers,4–8 the predictive value of transcriptomic signatures for Ipi/Nivo and nivolumab monotherapies remains limited, as evidenced by trials such as IMmotion150 and JAVELIN Renal 101.5 9–11

Endogenous retroviruses (ERVs), which constitute 8% of the human genome, are predominantly silenced in somatic tissues,12 13 have been observed to be transcribed in various cancers.14 15 The reactivation of ERVs, linked to the induction of inflammatory responses, has been demonstrated to enhance ICI efficacy in preclinical models.15–17 Specifically, in ccRCC, the expression of certain ERVs has been correlated with tumor immune signatures and responsiveness to ICIs, although in a limited study of 24 patients.18

Our study aims to enhance the robustness of RNA-based biomarkers by performing extensive ERV profiling across a combined cohort of 229 patients. This analysis includes data from the CheckMate clinical trials, complemented by insights from the BIONIKK trial, which pioneers biomarker-driven patient selection based on the Descartes classifications.19 20 We propose an ERV stratification system that identifies two specific ERVs with divergent prognostic values. This dual ERV system seeks to propel the discovery of potent biomarkers and significantly improve patient outcomes in the realm of immunotherapy.

Materials and methodsRNA sequencing profiling

RNA was extracted from frozen samples employing the well-established AllPrep DNA/RNA kit (Qiagen, catalog# 80204) following the manufacturer’s instructions. The library preparation and sequencing processes adhered to rigorous protocols previously documented, which included total rRNA depletion and the use of a HiSeq4000 sequencer for paired-end 2×100 bp sequencing.21 Total RNA-seq is known to capture a more comprehensive array of ERV transcripts, including those often missed by polyadenylated RNA-seq approaches.22 23 Quality control of paired-end total RNA-seq FASTQ files from the BIONIKK trial was initially performed using FastQC, which assessed quality at the base and read levels. Subsequently, raw reads were preprocessed to remove adapters and low-quality sequences (Phred quality score below 20) using cutadapt, and reads shorter than 40 bases were discarded. Residual reads mapping to rRNA sequences were removed using bowtie2. Only samples meeting these quality standards were advanced for further analysis.21 Cleaned and filtered reads were then aligned to the GRCh38 reference genome using the STAR aligner to generate RNA-seq BAM files.24 Gene annotations were provided by GENCODE (Release 27). Gene-level quantification was achieved with STAR, computing fragments per kilobase million and transcripts per kilobase million metrics for expression analysis.

Human proviral endogenous retroviruses mapping

To quantify the expression of ERVs, we used ERVmap, a bioinformatics pipeline designed for the locus-specific quantification of ERV expression. This pipeline combines stringent filtering criteria for RNA-seq reads mapping to ERV loci with a reference annotation database. The database comprises a curated list of 3220 ERVs that resemble full-length proviral sequences, each with a unique chromosomal location identified as autonomous/proviral. These ERVs, primarily autonomous elements averaging 7.5 kb in length, represent a comprehensive resource for analyzing specific autonomous ERV genomic loci throughout the human genome.

For our analysis, paired-end total RNA-seq reads were aligned to these 3220 proviral loci in the human genome (hg38) using the Burrows-Wheeler Aligner.12 ERVmap applies a stringent filtering criterion to ensure high-fidelity mapping: each read must (1) match only one best locus, (2) have at least one more mismatch at the second-best match, and (3) be excluded if it has more than three mismatches. Following filtration based on default parameters, ERV-specific reads were quantified and normalized using DESeq2, and log2-transformed for subsequent analysis.25 To maintain focus on biologically significant features and to enhance statistical power, only ERVs with normalized counts greater than 1 in at least 10% of the samples were considered to have effective expression. Coordinates of all the unique elements were compared with a custom repeat region annotation previously built using RepeatMasker for GRCh38 (https://www.repeatmasker.org). Distinct regions of the same provirus were annotated separately using the splicejam R package.26

External transcriptomic profiles and ERV expression from CheckMate trials

Transcriptomic and ERV expression data for the CheckMate cohort were sourced from prior studies, involving 311 metastatic ccRCC tumors from patients participating in CheckMate-009, CheckMate-010, and CheckMate-025 trials.11 This dataset comprised 181 patients treated with nivolumab and 130 treated with everolimus, with survival data also retrieved. For the CheckMate data, employed as a secondary data source, ethical approval details and patient consent are as described in the original publication.11 We accessed the preprocessed ERV expression data for this cohort (DESeq2-normalized and log2-transformed), quantified against a reference ERV sequence set,27 of which the exact locus information was utilized by ERVmap.25 A total of 1717 ERVs were identified as expressed in at least one sample across the CheckMate cohorts independently. Of these, 1640 ERVs (96%) showed a high concordance rate with those identified by the ERVmap algorithm. Ultimately, a subset of 666 ERVs with substantial expression, intersecting between the BIONIKK and CheckMate cohorts, was selected for detailed analysis in this study.

Identification of prognostic ERVs through bootstrap approach

We adopted bootstrap methods as our primary analytical strategy to identify significant ERVs with robustness and reliability against data permutation. The process began by randomly selecting 90% of the samples without replacement from each cohort to generate bootstrap subsets. This subsampling was integral to our approach, allowing us to mitigate potential biases associated with specific subsets of data.

For each bootstrap subset, univariate Cox proportional hazards regression analyses were performed for each ERV. This procedure was repeated 1000 times for each cohort, which enabled us to estimate the distribution of HRs for each ERV and assess their stability across different permutations of the data. Throughout these iterations, we monitored how frequently each ERV demonstrated statistically significant HRs that were consistent, either indicating risk (HR>1) or protective effects (HR<1). An ERV was considered robust if its HR showed consistent stability in all iterations and statistical significance in over 700 of the 1000 iterations for both cohorts. This stringent threshold was set to ensure that only ERVs with highly reproducible effects were considered significant, thereby minimizing the risk of false positives due to random variation.

After establishing the robustness of the ERVs through bootstrap analysis, we recalculated the overall HR for each ERV using the full dataset from each cohort. This final step was crucial for confirming the predictive power and clinical relevance of the ERVs, as it leveraged the entire cohort to validate the findings obtained from the bootstrap subsets. This comprehensive approach ensures that the identified ERVs are not only statistically significant but also practically meaningful for predicting clinical outcomes in metastatic ccRCC.

DNA methylation profiling and quantification

Genomic DNA was extracted from frozen samples, as previously described.21 Following bisulfite conversion of genomic DNA, profiling of DNA methylation was done using the Illumina EPIC BeadArrays.21 We used the R package ChAMP to extract the raw signal intensities and calculate the β value for each probe.28 The β value was calculated as M/(M+U), where M and U, respectively, refer to the preprocessed mean methylated and unmethylated probe signal intensities.

Calculation of a DNA methylation-based epigenetic silencing index

We recently developed an index of DNA methylation-based epigenetic silencing (iMES) that is tightly associated with primary resistance to immunotherapy.29 Specially, we used iMES R package, with DNA methylation status as input (β value >0.2 as cut-off: 1 for methylated and 0 for unmethylated). Each iMES score is then normalized on a scale from 0 to 1 through a min-max normalization process, ensuring a standardized assessment across samples.

Estimation of TME cell abundance

The composition of the TME for each ccRCC case was estimated using the MCPcounter R package, which calculates abundance scores for specific immune and stromal cell populations.30 This included eight types of immune cells, including T cells, CD8+T cells, cytotoxic lymphocytes, natural killer (NK) cells, B cell lineage, monocytic lineage, myeloid dendritic cells, and neutrophils, and two stromal cell populations (endothelial cells and fibroblasts). To determine the functional orientation of TME, we used signatures derived from the literature,29 including T cell activation, T cell survival, regulatory T cells (Tregs), and myeloid-derived suppressor cell (MDSC). Scores for each signature were calculated as the geometric mean of signature expression.

Differential and enrichment analysis

Differential expression analysis used the MOVICS package with a limma-based approach. For gene set enrichment analysis (GSEA), a preranked gene list based on differential expressions was analyzed with clusterProfiler to identify hallmark pathway enrichments.31–33

Assessment of transcriptomic signatures

For samples with available RNA-seq profiles, scores were derived for three IMmotion150 signatures, including angiogenesis (Angio): VEGFA, KDR, ESM1, PECAM1, ANGPTL4, and CD34; T-effector (Teff): CD8A, EOMES, PRF1, IFNG, and CD274; myeloid inflammation (Myeloid): IL-6, CXCL1, CXCL2, CXCL3, CXCL8, and PTGS2.9 The 26-gene JAVELIN Renal 101 Immuno signature includes CD3G, CD3E, CD8B, THEMIS, TRAT1, GRAP2, CD247, CD2, CD96, PRF1, CD6, IL7R, ITK, GPR18, EOMES, SIT1, NLRC3, CD244, KLRD1, SH2D1A, CCL5, XCL2, CST7, GFI1, KCNA3, and PSTPIP1.5 The tumor inflammation signature (TIS) includes PSMB10, HLA-DQA1, HLA-DRB1, CMKLR1, HLA-E, NKG7, CD8A, CCL5, CXCL9, CD27, CXCR6, IDO1, STAT1, TIGIT, LAG3, CD274, PDCD1LG2, and CD276.34 Signature scores were calculated as the median value of z-scored expression for the constituent transcripts. For CheckMate cohort, precalculated scores of Angio, Teff, Myeloid, and Immuno were collected directly from the prior study.11

Statistical analysis

Comprehensive statistical analyses were executed in R V.4.2.2. The Mann-Whitney U test and Fisher’s exact test were employed for continuous and categorical data, respectively, with missing data excluded. Survival outcomes were assessed via Kaplan-Meier curves and differences were evaluated using the log-rank test, with post hoc pairwise comparisons adjusted for false discovery rate (FDR) using the Benjamini-Hochberg method. Cox proportional hazard regression was used to calculate HRs and 95% CIs. Model performance was assessed through time-dependent receiver operating characteristic (ROC) analysis using the survivalROC package, which determined the area under the curve (AUC) at various time points. The Youden Index was used to optimize cut-offs, balancing sensitivity and specificity to enhance model interpretation and predictive accuracy. Brier scores were also calculated to assess prediction accuracy, representing the average squared differences between observed survival statuses and predicted probabilities, with values ranging from 0 (perfect prediction) to 1. All unadjusted statistical tests were considered significant at a p<0.05.

ResultsPatient characteristics

Our analysis encompassed 229 patients with available RNA-seq data undergoing ICI therapy for metastatic ccRCC across two independent cohorts (figure 1A). Within the BIONIKK trial, RNA-seq data were available for 48 patients who received first-line Ipi/Nivo combination therapy. Similarly, from the CheckMate trials, 181 patients, predominantly treated with nivolumab monotherapy as second or subsequent lines of treatment, had RNA-seq data available for analysis. The median age was 62 years (IQR: 54–69), with a predominance of male patients (79%). International Metastatic RCC Database Consortium (IMDC) risk stratification classified 40 patients (24%) as favorable risk and 124 (76%) as intermediate or poor risk. Among 207 samples assessed for sarcomatoid differentiation, 25 (12%) exhibited this feature. Based on Response Evaluation Criteria in Solid Tumors (RECIST) V.1.1 criteria, 67 patients (29%) achieved a complete (CR) or partial response (PR), whereas 153 (67%) experienced stable (SD) or progressive disease (PD); 9 cases were not evaluable for response. We defined long-term ICI responders as those patients achieving a CR, PR, or SD with a progression-free survival (PFS) exceeding 6 months. Conversely, non-responders were characterized as patients experiencing PD, or those with SD or PR where progression occurred within 6 months. Additionally, for the CheckMate cohort, in cases where RECIST was not evaluated (NE), we empirically classified patients who died within 6 months as non-responders. Conversely, patients who were alive without documented progression at the 6-month mark were considered responders. Out of the 229 patients, 106 (46%) were identified as demonstrating a long-term response, while 123 (54%) did not.

Figure 1

Overview of clinical trial cohorts and ERV analyses. (A) Comparative clinical characteristics of ccRCC patients in the CheckMate and BIONIKK trials. (B) Illustration of the chromosomal distribution for 666 ERVs identified with consistent expression across both cohorts. (C) Bootstrap consensus dot plot displays the congruence of bootstrap results across both cohorts. Dots along the diagonal represent ERVs with similar bootstrap outcomes in both cohorts. The size of each dot corresponds to the frequency of an ERV being prognostically significant in 1000 bootstrapping iterations, while the color indicates the stability of the HR determined by univariate Cox regression: blue for protective association (HR<1 in all cases), red for risk association (HR>1 in all cases), and mixed colors for variability in prognostic relevance across data permutations. The three ERVs highlighted in the top right corner denote those with stable HR and consistent univariate significance in over 70% of bootstrap runs. (D) A cytoband plot delineates the chromosomal positions and prognostic significance of 666 ERVs, highlighting 3 ERVs with stable prognostic value across cohorts using full dataset as determined by univariate Cox regression. Blue lines represent ERVs associated with a protective effect (HR<1, p<0.05), and red lines denote ERVs linked to a higher risk (HR>1, p<0.05). ERVs marked with a green circle retain their statistical prognostic significance after refinement by secondary clinical endpoint. ERV, endogenous retrovirus.

Association between ERVs and ICI efficacy

A total of 666 ERVs exhibited effective expression across both cohorts (figure 1B). Bootstrap analysis using univariate Cox proportional hazard regression, with PFS as the outcome, robustly identified three prognostic ERVs across both cohorts. These included E4421_chr17, associated with increased PFS risk, and two protective ERVs, E1659_chr4 and E1921_chr5 (figure 1C). Subsequent validation using the full dataset confirmed the predictive significance of these ERVs, ensuring robustness beyond the bootstrap subsets (online supplemental tables S1–S2). This step leveraged the larger sample size to refine HR estimates, providing a stronger statistical foundation for our findings and affirming their applicability across the entire study population.

Further analysis of secondary outcomes—overall survival (OS) for CheckMate and time-to-second treatment (TST) survival for BIONIKK—refined these ERVs’ prognostic significance (figure 1D). Notably, E4421_chr17 was marginally associated with worse OS in CheckMate (p=0.094; online supplemental table S3) and significantly correlated with poorer TST in BIONIKK (p=0.043; online supplemental table S4). Among the protective ERVs, E1659_chr4 showed superior prognostic relevance for both OS (p=0.043) and TST (p=0.054) compared with E1921_chr5 (OS: p=0.116, TST: p=0.024) (online supplemental table S3–S4), leading to a refined dual ERV panel for predicting ICI efficacy in metastatic ccRCC.

Dual ERVs analysis refined four distinct metastatic ccRCC subgroups

Median expression level of dual ERVs (normalized and log-transformed count values of 7.33 for E4421_chr17 and 4.95 for E1659_chr4) stratified the CheckMate cohort into four subgroups (ERV4): H17L4, H17H4, L17L4, and L17H4. “H” denotes high expression and “L” low expression, with the numerals referencing the corresponding chromosomes—“17” for E4421 on chromosome 17 and “4” for E1659 on chromosome 4 (eg, H17L4 indicates high expression of E4421_chr17 and low expression of E1659_chr4). The L17L4 subgroup defined the largest group of cases (n=54; 30%), followed by H17L4 (n=47; 26%), H17H4 (n=44; 24%) and L17H4 (n=36; 20%) (figure 2A).

Figure 2

ERV classification and its association with immune checkpoint inhibitors outcomes. (A) Stratification of the CheckMate cohort into four ERV subgroups (ERV4) based on the expression of key ERVs (E4421_chr17 and E1659_chr4), correlated with long-term ICI response categories. (B) Kaplan-Meier survival curves for the ERV4 in the nivolumab arm of the CheckMate cohort, analyzed by the log-rank test. (C) Sankey plot depicting the relationship between ERV4 subgroups, three-tier ERV risk model (ERV3), RECIST outcomes, and long-term ICI response. Accompanying stacked bar plots show the distribution of long-term ICI responses within each ERV3 category in the CheckMate cohort. (D) Time-dependent ROC curves of the numerical representation form of ERV3 in CheckMate cohort over a 10-year period, indicating the true positive rate (TPR) versus the false positive rate (FPR). (E–G) Same as (A–C) but for Ipi/Nivo arm of the BIONIKK cohort. (H) Time-dependent ROC curves of ERV3 within 12 months post-treatment in BIONIKK cohort. (I) Forest plot illustrating HRs from univariate (above dashed line) and multivariate (below dashed line) analyses, highlighting significant prognostic factors (in red) with p<0.05. Each variable’s comparison group is indicated in parentheses, with HRs computed relative to the reference group. For instance, “Age (>60)” indicates that the HR is calculated for patients older than 60, compared with those 60 or younger. AUC, area under the curve; CR, complete response; ERV, endogenous retrovirus; FDR, false discovery rate; ICI, immune checkpoint inhibitor; IMDC, International Metastatic RCC Database Consortium; PD, progressive disease; PFS, progression-free survival; PR, partial response; ROC, receiver operating characteristic; SD, stable disease; TST, time-to-second treatment.

The analysis of clinical characteristics across ERV4 demonstrated significant associations with treatment outcomes (table 1), particularly showing strong correlations with both PFS (p=0.0005) and OS (p=0.0029) (figure 2B). The post hoc analysis highlighted the L17H4 subgroup as exhibiting the most favorable outcomes, significantly outperforming H17L4 in median PFS (5.9 vs 1.9 months, p=0.0003, FDR=0.002) and OS (39.1 vs 15.9 months, p=0.002, FDR=0.010) (figure 2B, online supplemental figure S1). The disparities between H17L4 and other subgroups (eg, H17H4 and L17L4) further emphasize the robustness of the ERV-based stratification in identifying patients at higher risk of poorer outcomes under nivolumab monotherapy (FDR<0.05) (online supplemental figure S1). Adjusted comparisons involving H17H4 revealed no significant differences with L17H4 or L17L4, indicating overlapping survival characteristics among these subgroups.

Table 1

Comparisons of clinical characteristics across ERV four subgroups in CheckMate and BIONIKK cohorts

As L17L4 and H17L4 showed generally intermediate survival, indicating an ordinal nature of the ERV4 with varying prognostic implications, we introduced a numeric transformation to a three-tier risk model (ERV3) for analytical clarity: H17L4 labeled as one for high risk, L17L4/H17H4 as two for intermediate risk, and L17H4 as three for low risk. This helped clarify the gradient in survival outcomes, which was further supported by doubled ICI long-term response rates in the low-risk L17H4 (50%) compared with the high-risk H17L4 (21%) (p=0.006) (figure 2C). Time-dependent ROC analysis underscored the discriminative power of the ERV risk model, with AUC values of 0.76 for PFS and 0.71 for OS at 10 years (figure 2D).

Applying the same median cut-off values from CheckMate to the BIONIKK cohort validated the ERV4 framework in a different clinical setting (table 1, figure 2E), highlighting significant distinctions in both PFS (p=0.012) and TST (p=0.046) for patients treated with Ipi/Nivo combination immunotherapy. Particularly, the L17H4 subgroup in BIONIKK experienced markedly extended median PFS and TST of 23.2 and 20.8 months, respectively, starkly contrasting with outcomes in the H17L4 subgroup, where median PFS was 4.2 months and TST was 13.6 months (figure 2F). Post hoc pairwise comparisons confirmed significant survival advantages for L17H4 over H17L4 (PFS: p=0.001, FDR=0.007; TST: p=0.003, FDR=0.018) (online supplemental figure S1). Conversely, comparisons involving L17H4 or H17L4 with other subgroups, or comparisons between other subgroups did not show significant differences post-adjustment.

Moreover, converting the ERV4 into ERV3 in BIONIKK revealed that nearly all patients (94%) in the low-risk L17H4 showed a long-term clinical response to ICI, again doubling the response rate compared with 43% in the high-risk H17L4 group (p=0.004, figure 2G). These generally higher long-term response rates in the BIONIKK cohort compared with those in the CheckMate cohort underscore the efficacy of Ipi/Nivo over nivolumab monotherapy and demonstrate the discriminative ability of ERV3 in predicting PFS and TST for patients treated with combination immunotherapy, with average AUC values of 0.74 and 0.71, respectively (figure 2H).

The Brier score analysis for the BIONIKK and Checkmate cohorts generally indicates low predictive error of ERV3 (online supplemental figure S2). In BIONIKK, both PFS and TST exhibited increasing Brier scores over time, indicating challenges in accurately predicting longer-term outcomes. Notably, ERV3 was particularly effective in predicting short-term outcomes for TST, with Brier scores remaining below 0.05 within the first 3 months post-treatment. Conversely, in the Checkmate cohort, while Brier scores for PFS remained consistently low, scores for OS increased with time and stabilized around 0.2 after 2 years, suggesting the model’s difficulty in accurately forecasting long-term mortality.

Noteworthy is the absence of statistical correlation between the ERV4 classifications and patient prognosis following treatment with sunitinib in the BIONIKK cohort or everolimus in the CheckMate cohort, suggesting specificity in the predictive value of ERV-based classification system for immunotherapy outcomes (all, p>0.05; online supplemental figure S3).

ERV risk model as an independent prognostic factor

We then assessed the independent prognostic value of ERV3 risk model, multivariate analyses were conducted to adjust for major clinical prognostic features across both cohorts. In the CheckMate cohort, after adjusting for other prognostic factors, ERV3 remained a significant predictor of PFS (p=0.0002) and OS (p=0.032) (figure 2I). Similarly, in the BIONIKK cohort, ERV3 was independently associated with PFS (p=0.024), demonstrating its utility as a prognostic marker. While the association with TST approached statistical significance (p=0.051) (figure 2I), this was potentially limited by the smaller sample size of the BIONIKK cohort. Despite this, the data reinforced the close association between ERV3 and ICI efficacy in patients with metastatic ccRCC. Moreover, this implies that ERV3 risk categories could enhance the IMDC risk score’s accuracy in forecasting outcomes, especially in the context of immunotherapy.

Association between ERV risk model and TME

To better elucidate the relationship between ERV3 and ICI responses, we examined the abundance of 10 cellular components, the functional orientation and the expression of immune checkpoint genes within the TME for both cohorts across the three risk categories with different prognosis (figure 3A,B).

Figure 3

ERV risk model and their relationship with the TME. (A) Box plots depict the distribution of ten TME cellular components across ERV3 model (low-risk L17H4, intermediate-risk L17L4 or H17H4, and high-risk H17L4) in the CheckMate cohort, with scaled mean abundances annotated as heatmap below. (B) Same as (A) but for BIONIKK cohort. (C, D) GSEA plot showing pathways of interest upregulated in the low-risk L17H4 compared with the high-risk H17L4, for both CheckMate and BIONIKK cohorts, respectively. *p<0.05 by Mann-Whitney U test. ERV, endogenous retrovirus; FDR, false discovery rate; TME, tumor microenvironment.

Although variations in TME composition among the ERV3 were noted, a consistent pattern emerged: the low-risk L17H4 was enriched with endothelial cells in comparison to the high-risk H17L4 (CheckMate: p=0.008; BIONIKK: p=0.033) (figure 3A,B). This finding was further substantiated through GSEA, which highlighted a significant activation of hypoxia hallmark pathways in L17H4 versus H17L4 across both cohorts (figure 3C,D). Additionally, in the CheckMate cohort specifically, the angiogenesis hallmark pathway was notably activated in L17H4 (figure 3C).

The abundance of immune cell infiltration presented variations between the two cohorts, with the exception of myeloid dendritic cells, which were observed in greater abundance in the high-risk H17L4 (CheckMate: p=0.033; BIONIKK: not significant). Intriguingly, in CheckMate cohort, immune cells seemed to be depleted in the low-risk L17H4, while the BIONIKK cohort displayed an increased infiltration of several immune cells, such as NK cells and neutrophils in the low-risk category (figure 3A,B). While these findings did not reach statistical significance—likely due to the limited sample size—they align with recent reports suggesting these factors as potential predictors of ICI response.35 36 Furthermore, the high-risk H17L4 may harbor a slightly greater presence of immunosuppressive agents such as Tregs and MDSCs, alongside elevated levels of exhaustion markers (CTLA4, PDCD1, LAG3) compared with the low-risk L17H4 in the CheckMate cohort (all, p<0.1), while statistical significance was not reached in BIONIKK cohort (online supplemental figure S4).

Moreover, we revealed that the inflammatory response was activated in both cohorts, and activation of interferon-γ hallmark pathway was specifically observed in the BIONIKK cohort (figure 3C,D). This underscores the intricate and multifaceted interactions among ERV3 risk categories, cellular composition of the TME, and the consequent functional responses to ICIs.

ERV risk model surpasses transcriptomic signatures in predicting ICI efficacy

In our investigation into the prognostic value of established transcriptomic signatures for ICI response, we assessed four widely recognized signatures within the CheckMate and BIONIKK cohorts: Teff, Myeloid, Immuno, and TIS.5 9 34 Notably, within the CheckMate cohort, these signatures did not show a statistically significant correlation with clinical outcomes in patients undergoing nivolumab monotherapy (all, p>0.05; figure 4A). As vasculature might influence the ICI efficacy in ccRCC, evidenced by GSEA, we also examined the Angio signature.9 29 Unlike the others, the Angio signature demonstrated a significant association with improved OS (p=0.002), but not with PFS (p=0.553). Further multivariate analysis of CheckMate underscored the ERV3’s independence as a prognostic factor for enhanced outcomes, even when accounting for the Angio signature’s influence on OS (HR=0.80, 95% CI: 0.69 to 0.92, p=0.002). Of note, among patients receiving the Ipi/Nivo combination therapy in BIONIKK, only ERV3 emerged with prognostic significance, predicting both PFS (HR=0.50, 95% CI: 0.31 to 0.81, p=0.005) and TST (HR=0.48, 95% CI: 0.28 to 0.82, p=0.008) in univariate analyses. This finding was exclusive, as no other transcriptomic signatures showed statistical relevance to patient outcomes (all, p>0.05; figure 4A).

Figure 4

Comparative analysis of ERV risk model and transcriptomic signatures in predicting immunotherapy response. (A) Forest plot shows HRs from univariate (above dashed line) and multivariate (below dashed line) analyses, with significant predictors (in red) indicated by p<0.05, highlighting the predictive value of ERV3 relative to established transcriptomic signatures. (B) Time-dependent AUC, specificity, and sensitivity for ERV3 and other signatures, evaluated at the Youden index optimal cut-off in the CheckMate cohort. Performance averages are calculated across selected time points. (C) Same as (B) but for the BIONIKK cohort. AUC, area under the curve; OS, overall survival; PFS, progression-free survival; TIS, tumor inflammation signature; TST, time-to-second treatment.

In our comprehensive evaluation focusing on the predictive performance using time-dependent AUC, sensitivity and specificity at the Youden index optimal cut-off, the ERV3 risk model demonstrated superior predictive power compared with other signatures across various time points. Particularly through AUC analysis, ERV3 consistently showed higher values, suggesting a robust capability to distinguish between patient outcomes more accurately than other markers (figure 4B,C). However, a more nuanced analysis of sensitivity and specificity revealed patterns that reflect the complex interplay between ERV expressions, and the types of immunotherapies administered. In the CheckMate cohort, where patients received nivolumab monotherapy, ERV3 exhibited higher specificity (figure 4B). This suggests that ERV3 is particularly effective in accurately identifying patients who are less likely to benefit from nivolumab, thereby avoiding unnecessary treatment. Conversely, the lower sensitivity observed indicates a limitation in identifying all patients who would benefit from the therapy. In contrast, within the BIONIKK cohort, where patients were treated with the Ipi/Nivo combination immunotherapy, ERV3 showed a tendency toward higher sensitivity (figure 4C). This suggested a stronger ability of ERV3 to identify responders to combination therapy, capturing a broader spectrum of patients likely to benefit. However, the specificity observed was generally lower, particularly when considering TST, suggesting that the ERV3 was less precise in determining who would not require subsequent treatment imminently.

Refinement of immunotherapy response prediction via the integrated epigenetic-retroviral score

Emerging evidence underscores a complex interplay between ERVs and epigenetic alterations, highlighting their significant roles in modulating immunotherapy responses.37 38 Building on our prior discovery of an epigenetic silencing marker, iMES, linked to immunotherapy resistance in ccRCC,29 we aimed to enhance the predictive accuracy of immunotherapy outcomes by integrating insights from ERV3 risk model with epigenetic markers.

In the analysis of 48 patients from the BIONIKK cohort treated with Ipi/Nivo, we evaluated DNA methylation profiles in 27 cases. Univariate analysis indicated a marginal association of iMES with resistance to immunotherapy, whereas ERV3 showed a significant correlation with the response to treatment (table 2). Notably, on adjusting for ERV3, the association of iMES with immunotherapy resistance became significantly more pronounced (table 2). This observation suggests that ERV3 may act as a confounding factor, influencing the relationship between iMES and immunotherapy outcomes. This enhancement led us to propose the integration of ERV3 with iMES, hypothesizing that such a combined approach could offer superior predictive insights into immunotherapy efficacy.

Table 2

Impact of ERV three-tier risk model (ERV3) and iMES on survival outcomes: univariate and multivariate Cox regression analyses in the BIONIKK cohort

To operationalize this hypothesis, we introduced the Integrated Epigenetic-Retroviral (IER) score, a prognostic index calculated by dividing the iMES value by the numerical representation of the ERV3 risk model. This formula encapsulates both epigenetic and retroviral dimensions into a singular metric. The IER score significantly correlated with increased risk of immunotherapy resistance (figure 5A). Of note, this methodological integration significantly improved the association with immunotherapy resistance as time-dependent ROC analysis corroborated the IER’s superior predictive performance over individual ERV3 or iMES assessments within 6–18 months post-treatment (all, p<0.05; figure 5B).

Figure 5

Enhancement of immunotherapy response prediction through an IER score. (A) Kaplan-Meier survival analysis comparing patients with high versus low IER scores in the BIONIKK cohort, evaluated using the log-rank test with HRs derived from univariate Cox regression. (B) Time-dependent AUC demonstrates the superior predictive accuracy of the IER score than using either the ERV3 or the iMES epigenetic marker alone, assessed from 6 to 18 months following immunotherapy treatment in the BIONIKK cohort. *p<0.05, **p<0.01 by Mann-Whitney U test. IER, Integrated Epigenetic-Retroviral; iMES, index of DNA methylation-based epigenetic silencing; mAUC, mean area under the curve; PFS, progression-free survival; TST, time-to-second treatment.

Discussion

Human ERVs, stemming from ancient viral incorporations into the human genome, are typically dormant but can become reactivated in cancer cells, often due to the hypomethylation which is common in oncogenesis.39 40 This re-expression, while contributing to genome instability and potentially driving oncogenesis,41–43 also plays a dual role in the immune landscape of cancer—on one hand, producing tumor-associated neoantigens that may stimulate an antitumor immune response, and on the other, potentially exerting immunosuppressive effects through mechanisms that include the modulation of immune cell activation.44–46 Our study advances the understanding of ERVs as predictive biomarkers for immunotherapy in metastatic ccRCC, particularly emphasizing the significance of ERVs like E4421_chr17 and E1659_chr4. These ERVs may enhance or suppress the TME response to ICIs, which align with emerging evidence that underscores the role of viral mimicry in modulating the TME and influencing the efficacy of immunotherapy—an insight that conventional transcriptomic signatures fail to capture.47

E4421_chr17 is classified as containing elements from both the long-interspersed element (LINE) and long terminal repeat (LTR) classes, specifically from the ERVK and L1 families. HERV-K expression has been frequently documented to increase in tumor cells. Importantly, HERV-K-derived proteins (eg, Rec and Np9) can act as tumor-specific antigens, a class of neoantigens, and induce immune responses in different types of cancer.48 In ccRCC, the expression of ERVs like E4421_chr17 from the ERVK family correlates with increased immune checkpoint activation and altered immune cell dynamics, potentially leading to an environment conducive to immune evasion. Chronic exposure to these antigens might contribute to immune tolerance or T-cell exhaustion, with increased recruitment of Tregs and MDSCs, as evidenced in our findings. Moreover, ERVs modulate the immune response by interacting with cytosolic nucleic acid sensors such as RIG-I and STING, enhancing inflammation and potentially leading to fibroinflammation in diseases like kidney disease, illustrating the broader implications of ERV activity beyond oncogenesis.49

E1659_chr4, also containing elements from the LINE and LTR classes, is associated with the ERV1, ERVL-MaLR, and L1 families. Previous studies reported several ERV loci from ERV1 family associated with cytotoxic T cell presence and antigenic epitopes and response to ICI in ccRCC,15 18 while none of them were identified as prognostic in first-line immunotherapy. Recent studies in other cancers have shown that antibodies against ERV proteins can enhance the effectiveness of immunotherapy by promoting tumor immunity.50 For instance, in lung adenocarcinoma, local B cell responses against ERV envelope glycoproteins were amplified by ICI, enhancing antitumor activity.37 These observations align with our findings that ERV expression can modulate immune responses in ccRCC, suggesting broader implications of ERVs in cancer immunotherapy across different tumor types.

Our findings advocate for an ERV-based classification as a potent prognostic tool, standing independently from other clinical factors and potentially augmenting the IMDC risk scores that may not fully reflect the varied responses to treatment.51 Insights from examining the TME’s relationship with ERV risk model elucidate how differential responses to ICI arise, with a pattern of endothelial cell presence emerging as a common denominator in successful treatments.29 Such insights align with observations from the BIONIKK trial, which reported higher response rates to Ipi/Nivo treatment in the proangiogenic group compared with those receiving sunitinib—a finding echoed by the IMmotion 150 trial’s identification of an angiogenesis signature predictive of response to combined atezolizumab and bevacizumab.5 9 20 Additionally, a study in melanoma linked greater vascular density with better outcomes from anti-CTLA-4 therapy.52 Furthermore, the myeloid compartment of the TME is highly heterogeneous, consisting of macrophages, dendritic cells, and MDSCs, each playing roles that can either support or inhibit tumor growth and immune evasion.53 This diversity and plasticity make myeloid cells crucial in determining the response and resistance to immunotherapy.

Crucially, ERVs, typically silenced by promoter DNA methylation, can become activated across various cancer types.54 Our study introduces the IER score, which combines ERV risk model with epigenetic silencing markers to enhance the prediction of immunotherapy responses. This innovative metric underscores the intricate relationship between genomic and epigenetic landscapes in modulating treatment efficacy, a concept that is beginning to gain traction in oncology research. By leveraging RNA and DNA from standard biopsies, our dual-layer analysis—combining ERV and methylation data—fits within current clinical workflows, enhancing predictive accuracy while remaining practical for clinical use. This approach marks a significant step toward comprehensive precision oncology, integrating multiple molecular dimensions to optimize patient-specific treatment strategies.

The limitations of this study should be acknowledged alongside its strengths. Although the analyses conducted were post hoc and not predesigned, we used reliable data from two specific prospective trials, reducing the risk of collection bias often associated with retrospective data. However, there is potential selection bias due to the reliance on analyzable tumor tissue from a highly selected cohort, which emphasizes the need for comparisons with the overall cohorts to enhance generalizability. Additionally, the variability in lines of therapy among patients with metastatic ccRCC could influence observed