Endogenous retrovirus group FRD member 1 is a potential biomarker for prognosis and immunotherapy for kidney renal clear cell carcinoma

Abstract

Introduction:

The activation of endogenous retroviral (ERV) genes in kidney renal clear cell carcinoma (KIRC) suggests the necessity for further research on their functions.

Methods:

In this study, KIRC and healthy cohorts were obtained from TGGA and GEO datasets. Subsequently, differential analysis and functional annotation were conducted using GO, KEGG, and GSEA. Clinical outcomes were then observed and utilized in the development of a nomogram.

Results:

We observed the general low expression of ERVFRD-1 in KIRC tumors compared to normal tissue (P < 0.001) across multiple cohorts. Differential analysis and functional annotation using GO, KEGG, GSEA analysis revealed significant involvement of ERVFRD-1 in tumor immunoregulation: a close relation to the infiltration levels of mast cells and Treg cell (P < 0.001) and occurrence with a variety of immune markers. Methylation status was then applied to uncover potential mechanisms of ERVFRD-1 in KIRC. Notably, higher expression levels of ERVFRD-1 were associated with extended overall survival, disease-specific survival, and progression-free survival. Finally, based on Cox regression analysis, we constructed a nomogram incorporating ERVFRD-1, pathologic T, and age, which exhibited promising predictive power in assessing the survival outcomes of KIRC patients.

Discussion:

To sum up, our study suggests that ERVFRD-1 plays a role in regulating immunological activity within the tumor microenvironment and is associated with overall survival in KIRC patients. ERVFRD-1 may therefore be a sensitive biomarker for diagnosis, immunotherapy, and prognosis assessment of KIRC.

1 Introduction

According to the latest GLOBOCAN report, kidney cancer is ranked 16th among the 36 types of cancers worldwide, with an estimation of over 431,000 new cases (accounting for 2.2% of all cancer types) and approximately 179,000 deaths (1.8% of all cancer-related deaths) (Sung et al., 2021). Kidney cancer is typically divided into three histological types: clear cell (TCGA-KIRC, WHO-ccRCC), papillary (TCGA-KIRP, WHO-ChRCC), and chromophobe (TCGA-KICH, WHO-PRCC) (Moch et al., 2022). KIRC is the most prevalent type, representing over 80% of all kidney cancer cases. Unfortunately, KIRC is highly malignant and not responsive to most chemotherapy or radiotherapy (Makhov et al., 2018). While surgical resection is currently an effective treatment for early-stage localized KIRC (Porta et al., 2019), around 25% to 50% of primary patients experience recurrence within five years following nephrectomy, and one-third of patients develop metastases (Kotecha et al., 2019). Therefore, it is urgent to discover novel and valuable biomarkers for early diagnosis and treatment of this particular type of kidney cancer.

Human endogenous retroviruses (HERVs) are estimated to make up approximately 8% of the human genome (Geis and Goff, 2020). These retroviruses are typically inherited through the germ line. Due to multiple mutations in their genomes, HERVs are usually defective and rarely capable of producing viral particles. Despite this, HERVs-derived proteins can significantly influence the regulation of gene expression and cellular functions in the host. Apart from viral transcripts, ERV-encoded proteins translated from intact open reading frames (ORFs) have been implicated in both pathological and physiological processes (Turner et al., 2001). Numerous studies have demonstrated that HERV families are transcribed in human tissues during various stages of development, in both healthy and diseased conditions (Suntsova et al., 2015). HERVs have therefore been associated with several biological processes, including placentation, stemness, aging, immune responses, neurodegeneration, autoimmunity, and carcinogenesis (Okahara et al., 2004; Balada et al., 2009; Zhang et al., 2019; Nair et al., 2021; Rangel et al., 2022).

One particular HERV, ERVFRD-1, produces a protein called Syncytin-2, which plays a significant role during embryogenesis by contributing to placental formation and facilitating immune tolerance of the maternal system towards the developing fetus. In an in vitro experiment, it was observed that MCA205 cells transfected with syncytin-2 exhibited long-lasting tumorigenic potential when engrafted into mice (Mangeney et al., 2007). Moreover, increased expression of ERVFRD-1 in cancers from the Cancer Genome Atlas (TCGA) with ERVFRD-1 fusions suggested that ERVFRD-1 promoted tumor growth by inhibiting the anti-tumor immune response of the host mice. This finding suggests that enhanced expression of ERVFRD-1 in human cancer cells could potentially contribute to tumor growth as well (Namba et al., 2021). However, the precise role of ERVFRD-1 in KIRC has not been investigated.

Cancer immunotherapy has emerged as a potential alternative or complement to traditional cytotoxic chemotherapy and radiotherapy. However, the effectiveness of immune checkpoint blockades (ICBs), a key approach in immunotherapy, is significantly influenced by the tumor microenvironment (Li et al., 2021). KIRC patients often exhibit higher immune-related scores. Interestingly, studies have linked HERVs expression to response to ICBs in certain cancer types. In both human and mouse lung adenocarcinomas, researchers have discovered that tumor-binding antibodies and local germinal center responses are elicited, with ERV-envelope glycoproteins identified as the primary anti-tumor antibody target. Additionally, ICBs are capable of amplifying the B cell responses within this process by targeting HERVs in both humans and mice (Ng et al., 2023). A different HERV pattern of association with the outcome of KIRC PD-L1 immunotherapy was observed by Lewis Au et al (Au et al., 2021). However, despite the known role of the ERVFRD-1 protein in suppressing immunity during placental development, its immunological impact in KIRC and its potential relationship with the effectiveness of immunotherapy remain unexplored. Therefore, further research is necessary to clarify the relationship between ERVFRD-1 and immunity in KIRC.

The objective of this study is to utilize bioinformatics to examine ERVFRD-1 expression and its clinicopathological and prognostic significance, as well as the underlying molecular mechanisms and immune cell infiltration in KIRC. By investigating these aspects, our study seeks to provide valuable insights that could assist clinicians in refining treatment strategies and enhancing outcomes for patients diagnosed with KIRC.

2 Materials and methods

2.1 Data collection and processing

We collected mRNA expression profiles and corresponding clinical characteristics of KIRC patients from two databases: the Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov) database (n = 532), and the Gene Expression Omnibus database (GEO, http://www.ncbi.nlm.nih.gov/geo/) (n = 144). The data obtained from these databases were in level 3 HTSeq-FPKM format and were normalized as transcripts per million reads (TPM). For a broader analysis encompassing multiple cancer types, we also retrieved RNA-sequencing data in TPM format from the University of California, Santa Cruz Xena (UCSC Xena website, https://xenabrowser.net/datapages/). The RNA-seq data of 33 pan-cancer samples including adrenocortical carcinoma (ACC), bladder cancer (BLCA), breast carcinoma (BRCA), cervical cancer (CESC), cholangiocarcinoma (CHOL), colorectal cancer (COAD), diffused large B-cell lymphoma (DLBC), esophageal cancer (ESCA), glioblastoma (GBM), head and neck squamous cell carcinoma (HNSC), kidney chromophobe (KICH), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), Acute Myeloid Leukemia (LAML), low-grade glioma (LGG), liver hepatocellular carcinoma(LIHC), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), mesothelioma (MESO), ovarian serous cystadenocarcinoma (OV), pancreatic adenocarcinoma (PAAD), pheochromocytoma and paraganglioma (PCPG), prostate cancer (PRAD), rectal cancer (READ), sarcoma (SARC), skin cutaneous melanoma (SKCM), stomach adenocarcinoma (STAD), testicular cancer (TGCT), thyroid carcinoma (THCA), thymoma (THYM), uterine corpus endometrial carcinoma (UCEC), uterine carcinosarcoma (UCS) and ocular melanomas (UVM), were obtained from TCGA. Normal samples were selected from the Genotype-Tissue Expression (http://commonfund.nih.gov/GTEx/) dataset.

2.2 Analysis of Differentially Expressed Genes (DEGs)

The cohort of patients with KIRC from TCGA was categorized into high and low ERVFRD-1 expression groups using the median ERVFRD-1 expression score. DEG analysis was performed between these two groups using the R package DESeq2 (Robinson et al., 2010; Love et al., 2014). The significance threshold was set at an adjusted P value of < 0.05, and a |log2-fold-change (FC)| > 1. We used Spearman’s correlation analysis to evaluate the association between ERVFRD-1 expression and the top 10 DEGs.

2.3 Protein-Protein Interaction (PPI) network analysis

Following the differential gene expression analysis, we constructed a PPI network using the Retrieval of Interacting Genes online database (STRING, https://www.string-db.org/). A confidence score threshold of > 0.7 was set, while other parameters were kept at their default values. The resulting PPI network was then visualized using Cytoscape software (version 3.9.1) (Lotia et al., 2013). Next, we used CytoHubba, a plugin within Cytoscape software, to identify the top 10 hub genes among these DEGs (Chin et al., 2014).

2.4 Functional enrichment analysis

To gain insights into the functional implications of the DEGs, we performed functional enrichment analyses, specifically Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. These analyses were conducted using the R package ggplot2 (version 3.3.6). Additionally, we employed gene set enrichment analysis (GSEA) with the R package clusterProfiler (version 4.2.1) (Subramanian et al., 2005; Yu et al., 2012). 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.

2.5 Immune infiltration analysis

To assess the extent of immune infiltration in KIRC, we utilized the “estimate” package in R to calculate the stromal, immune scores and ESTIMATE scores (Yoshihara et al., 2013), and conduct a comprehensive analysis involving 24 immune cell types, visualized by ggplot2[3.3.6]. The single-sample Gene Set Enrichment Analysis (ssGSEA) technique was then employed to evaluate the relative enrichment score of these immune cells. This analysis was performed using the GSVA R package (Bindea et al., 2013; Hänzelmann et al., 2013). Furthermore, we examined the correlation between the expression of ERVFRD-1 and these immune cells using Spearman’s correlation analysis. The Wilcoxon rank sum test was utilized to determine differences in immune infiltration levels between groups with high and low ERVFRD-1 expression. Additionally, to examine the interaction between tumor immunity and markers, we utilized the ggstatsplot R package to analyze the correlation between the enrichment score of tumor mutation burden (TMB) and ERVFRD-1 expression. Then potential ICB response was predicted using the TIDE algorithm (Jiang et al., 2018), visualized by ggplot2(v3.3.3) and ggpubr(0.4.0).

2.6 DNA methylation analysis

To explore the underlying mechanism of ERVFRD-1 in KIRC, we employed the MethSurv database (https://biit.cs.ut.ee/methsurv/) to evaluate the DNA methylation status of ERVFRD-1 and prognostic value of ERVFRD-1 methylation levels. MethSurv is an online tool that allows for multivariable survival analysis based on DNA methylation data (Modhukur et al., 2018).

2.7 Survival analysis

We performed survival analysis using the Kaplan-Meier method and assessed statistical significance with the log-rank test, where the median expression level of ERVFRD-1 was used as the cut-off value. To evaluate the impact of clinical variables on patient outcomes, we conducted univariate and multivariate Cox regression analyses. In the univariate Cox regression analysis, prognostic variables with a significance level of P < 0.1 were included in the multivariate Cox regression analysis. The results were visualized in a forest plot using the R package ggplot2.

2.8 Construction and validation of the nomogram

We created a nomogram using the results from multivariate Cox analysis that included independent prognostic factors such as clinical characteristics and risk scores of KIRC patients. This nomogram aimed to predict the overall survival probability over 1-, 3-, and 5-year periods. The performance of the nomogram was evaluated using calibration plots while the concordance index (C-index) was utilized to quantify its discriminatory ability. The R package RMS (version 6.2-0) was utilized to create the nomogram and calibration plots. Additionally, we evaluated the predictive accuracy of the model using time-dependent receiver-operating characteristic (ROC) curves generated through the time ROC package.

2.9 Statistical analysis

R software (version 3.6.3) was used to conduct all statistical analyses. We evaluated the statistical significance of ERVFRD-1 expression in non-paired and paired tissues using the Wilcoxon rank-sum test and paired-sample t-test, respectively. The correlation between clinical features and ERVFRD-1 expression was assessed using both the Wilcoxon rank-sum test and logistic regression. All tests conducted were two-sided, and we regarded P values < 0.05 as statistically significant.

3 Results

3.1 Patient characteristics

We collected clinical and RNA-sequencing data for 532 patients from TCGA datasets, which consisted of 72 patients with matched normal tissue samples. The Supplementary Table S1 shows the clinicopathological characteristics of these patients. To augment the sample size, we acquired gene expression data from normal kidney tissues (n = 28) in the GTEx datasets. Furthermore, another dataset from the GEO dataset GSE53757 was introduced, which contained normal kidney tissues (n=71) and KIRC tissues (n=73).

3.2 Low expression of ERVFRD-1 associated with prognosis value in KIRC

Based on the pan-cancer analysis, ERVFRD-1 expression was generally low in various types of cancers such as BLCA, BRCA, CESC, CHOL, COAD, ESCA, KIRC, KIRP, LGG, LIHC, LUSC, PCPG, PRAD, READ, STAD, THCA, UCEC, and UCS (Figure 1A). Additionally, ERVFRD-1 expression was found to be significantly associated with prognosis in KIRC (Figure 1B). To further investigate the role of ERVFRD-1 in KIRC, we analyzed its expression in TCGA and GEO databases. As illustrated in Figure 1C, ERVFRD-1 expression was markedly lower in KIRC than in normal kidney tissues (p < 0.001). This trend was also observed in 72 paired kidney tissues (p < 0.001) (Figure 1D). Similarly, a significant decrease in ERVFRD-1 expression was found in KIRC compared to normal tissues in the GEO dataset GSE53757 (Figure 1E). Moreover, ERVFRD-1 expression demonstrated better predictive ability in distinguishing KIRC from normal tissues, exhibiting an area under the curve (AUC) of 0.952 (95% confidence interval [CI]=0.932-0.972) as indicated by the ROC curve (Figure 1F).

Figure 1

3.3 Association of ERVFRD-1 expression levels with clinical characteristics of KIRC

We utilized the TCGA database to investigate the impact of various pathologic features on ERVFRD-1 transcription in KIRC samples (Table 1). As displayed in Figure 2, ERVFRD-1 expression remained consistently lower in KIRC patients than in controls across different sub-groups. The sub-group analysis, which included T stage, N stage, M stage, pathologic stage, gender, race, age, histology grade, serum_calcium, hemoglobin, laterality, OS event, DSS event, and PFI event, also indicated a significant reduction in ERVFRD-1 expression in KIRC patients compared to controls. Moreover, the univariate logistic regression analyses showed statistically significant differences between the high and low ERVFRD-1 expression groups in terms of T stage (T3 & T4 vs. T1 & T2, OR = 0.414, 95% CI = 0.288 - 0.596, P < 0.001), M stage (M1 vs. M0, OR = 0.411, 95% CI = 0.245 - 0.688, P < 0.001), pathologic stage (Stage III & IV vs. Stage I & II, OR = 0.389, 95% CI = 0.272 - 0.557, P < 0.001), gender (Male vs. Female, OR = 0.604, 95% CI = 0.422 - 0.864, P = 0.006), race (white vs. Asian & Black or African American, OR = 0.368, 95% CI = 0.209 - 0.647, P < 0.001), and histologic grade (G3 & G4 vs. G1 & G2, OR = 0.449, 95% CI = 0.317 - 0.635, P < 0.001), respectively (Table 2).

Figure 2

3.4 Identification of DEGs in KIRC and PPI network analysis

Through a comparison of gene expression profiles between KIRC samples with high and low ERVFRD-1 levels, we discovered a total of 1348 DEGs, which consisted of 128 (9.50%) upregulated genes and 1220 (90.50%) downregulated genes (adjusted p-value < 0.05, |Log2-FC| > 1) (Figure 3A). We further explored the relationship between the top 10 DEGs (including AC090578.1, PAEP, FDCSP, SAA1, GOLGA6L7, APOA4, IGFBP1, RTL1, HHATL, and NMRK2) and ERVFRD-1, as illustrated in Figure 3B and Supplementary Table S2. We employed the online STRING tool to construct a PPI network that aimed to investigate the potential interactions within the DEGs (Supplementary Figure S1). Moreover, we identified a group of hub genes from the PPI network using the same tool. The network displayed a high level of complexity, and the top 10 hub genes were APOB, AHSG, FGA, FGG, FGB, APOC3, APOA4, IVL, SPRR1B, and RPTN, as shown in Figure 3C.

Figure 3

3.5 Functional enrichment analysis including GO, KEGG, and GSEA analysis

We then conducted GO and KEGG enrichment analyses to gain insights into the biological functions of the abovementioned DEGs. GO analysis revealed that these DEGs were primarily involved in “humoral immune response” and “complement activation, classical pathway” in terms of biological processes (GO-BP). In terms of cellular components (GO-CC), the DEGs were enriched in “immunoglobulin complex” and “immunoglobulin complex.” The molecular function (MF) analysis showed enrichment in “antigen binding” and “immunoglobulin receptor binding.” Additionally, the KEGG pathway analysis demonstrated enrichment in “neuroactive ligand-receptor interaction” and “complement and coagulation cascades.” These findings are depicted in Figure 3D and Supplementary Table S3.

Furthermore, we performed GSEA using GSEA/Molecular Signatures Database (MSigDB), which revealed significant enrichment in immune-related biological processes, including humoral immune response, immunoglobulin complex, and immunoglobulin receptor binding, among others. These results suggest that ERVFRD-1 may contribute to an enhanced immune phenotype in KIRC (Figure 4; Supplementary Table S4).

Figure 4

3.6 Correlation between ERVFRD-1 expression and immune infiltration in KIRC

To investigate the correlation between ERVFRD-1 expression and immune infiltrates in KIRC, we utilized the ESTIMATE algorithm to determine stromal scores, immune scores, and ESTIMATE scores. Our findings indicate that ERVFRD-1 expression is significantly and positively correlated with stromal score, but negatively correlated with immune score, suggesting a notable impact on stromal and immune cell infiltration (as shown in Figure 5A). Then we utilized ssGSEA and found a correlation between ERVFRD-1 expression and the levels of immune infiltrating cells, specifically mast cells (r=0.377, P < 0.001) and Regulatory T cells (Treg, r = -0.253, P < 0.001) (Figure 5B). The enrichment scores for mast cells were significantly higher in the group with high expression of ERVFRD-1 compared to the group with low expression, while Treg cells exhibited the opposite trend (all P < 0.001), as depicted in Figures 5C, D.

Figure 5

Additionally, we investigated the correlation between ERVFRD-1 expression and immune biomarkers, including immune modulators and tumor mutational burden (TMB). Our analysis revealed that the expression of ERVFRD-1 demonstrated a significant positive correlation with immunoinhibitors and most immunostimulators (e.g., CD70 and TNFSF14) (Figures 6A, B). In terms of TMB, higher TMB scores were observed in the low-expression group of ERVFRD-1, indicating the lower level of ERVFRD-1 was more likely to benefit from immunotherapy (P = 0.003) (Figure 6C). Additionally, a series of scattered plots, as depicted in Figures 6D, E, displayed the correlation between ERVFRD-1 with a few immunoinhibitors (ADORA2A, PD-L1, DIO1, KDR, KIR2DL1, KIR2DL3, PDCD1LG2, and TGFBR1) and immunostimlators (CD70, CXCL12, ENTPD1, IL6R, NT5E, RAET1E, TNFRSF18, and TNFRSF14). Furthermore, we examined 22 KIRC patients who had received immunotherapy (PD-L1/PD-1) from TCGA database and analyzed their response based on ERVFRD-1 expression levels. However, our results suggest no significant differences between high and low ERVFRD-1 expressing groups (Supplementary Figure S2).

Figure 6

3.7 Correlation between methylation and expression of ERVFRD-1 in KIRC

To shed light on the underlying mechanisms of ERVFRD-1 overexpression in KIRC tissues, we utilized the online tool “MethSurv database” to investigate the correlation between ERVFRD-1 expression levels and methylation status. Our analysis revealed that the DNA sequences of ERVFRD-1 contain two single Cytosine-phosphate-Guanine (CpG) sites, cg00966482 and cg26383454 (Figure 7A). Both of these methylation sites were observed to be methylated in KIRC, especially cg00966482. Furthermore, we identified patients who exhibited high levels of ERVFRD-1 methylation in the cg26383454 site had poorer overall survival rates when compared to those with lower levels of ERVFRD-1 methylation, while the methylation status of the cg00966482 site did not show any significant difference in survival between the high and low groups (Figures 7B, C).

Figure 7

3.8 Prognostic value of ERVFRD-1 in KIRC

We conducted survival analysis using the Kaplan-Meier method to evaluate the association between ERVFRD-1 expression and the survival outcomes of patients with KIRC. Patients were divided into high and low expression groups based on the median value of ERVFRD-1 expression. The high ERVFRD-1 expression group displayed a significantly better prognosis for overall survival (OS), disease-specific survival (DSS), and progression-free interval (PFI) in comparison to the low expression group (OS: hazard ratio [HR] = 0.59, 95% CI = 0.44 – 0.80, P < 0.001; DSS: HR = 0.37, 95% CI = 0.25 – 0.56, P < 0.001; PFI: HR = 0.48, 95% CI = 0.35 – 0.67, P < 0.001), as demonstrated by the corresponding Kaplan-Meier curves (Figures 8A–C).

Figure 8

Furthermore, we performed univariate and multivariate Cox regression analyses to identify prognostic factors. The multivariate analysis revealed that ERVFRD-1 expression (adjusted HR = 0.496, 95% CI = 0.293-0.839, P = 0.009), T2 stage (adjusted HR = 0.070, 95% CI = 0.005-0.935, P = 0.044), and age stage (adjusted HR = 2.050, 95% CI = 1.188-3.539, P = 0.010) were independent factors associated with OS in patients with KIRC (Figure 8D; Supplementary Table S5). Similarly, for DSS, ERVFRD-1 expression (adjusted HR = 0.375, 95% CI = 0.183 - 0.770, P = 0.008) was identified as a prognostic indicator (Supplementary Table S6). However, for PFI, ERVFRD-1 expression did not show prognostic significance (Supplementary Table S7).

We further assessed the prognostic value of ERVFRD-1 expression in different subgroups. Consistently, high expression of ERVFRD-1 was associated with favorable outcomes in various subgroups based on OS, including T3 & T4 stage, N0 stage, stages III & IV, grade 3& 4, gender, age > 60 years, white race, normal serum calcium, normal hemoglobin, and laterality (all P < 0.05) (Figure 9).

Figure 9

3.9 Construction and validation of a nomogram based on independent factors

We developed a nomogram using the independent prognostic factors for OS to provide a prognostic prediction tool for KIRC patients. The nomogram visually represented the total points assigned to each factor, where a higher total number of points indicated a worse prognosis (Figure 10A). Calibration curves were generated to assess the performance of the nomogram in predicting survival outcomes, demonstrating good agreement between predicted and observed survival probabilities (Figure 10B). The bootstrap corrected C-index of the nomogram was 0.710 (95% CI = 0.690 - 0.730), indicating a moderate predictive accuracy for OS in patients with KIRC. Additionally, we assessed the discriminative ability of ERVFRD-1 expression using time-dependent ROC curve analysis (Supplementary Figure S3). In summary, these observations suggested that the nomogram was an appropriate tool for predicting the prognosis of KIRC patients.

Figure 10

4 Discussion

Kidney renal cell carcinoma (KIRC), characterized by a microenvironment that is angiogenic, inflammatory, and highly immunogenic, has limited response to chemotherapy and radiation therapy (Kim et al., 2021). However, recent studies have shown the efficacy of immunotherapy as a treatment option for KIRC(Choueiri et al., 2020; Powles et al., 2022), underscoring the significance of investigating molecular markers that can provide prognostic information and insights into the immune response in this disease. The HERV family, remnants of ancient viruses, may not exhibit viral activity at the transcriptional level but has been found to enhance immune antigenicity. This characteristic has also been observed in tumor immunotherapy, making the exploration of HERVs crucial in understanding their role in the context of KIRC.

Evidence from previous publications has indicated abnormal expression of certain HERVs in urinary tract malignancies. For instance, upregulation of Syncytin-1, encoded by ERVW-1, is more frequently observed in bladder urothelial cancer compared to healthy tissues (Yu et al., 2014). A study by Veronika Weyerer and colleagues found significantly higher levels of expression for several HERVs, including Syncytin-1, ERV-H, ERV-Fc1, ERV-Fc2, ERV-T, ERV-Kenv, ERV-Kpol, ERV-E4-1, and ERV-E6q15 in KIRC lesions when compared to normal kidney samples (Weyerer et al., 2021). Our study, for the first time, provides direct evidence of a potentially vital role regarding the expression of ERVFRD-1 in KIRC. ERVFRD-1 is typically expressed in the placenta and thought to protect the fetus from rejection by the mother (Bolze et al., 2017). To begin with, our findings revealed from TCGA database that ERVFRD-1 exhibited low expression levels in KIRC tissues compared to normal tissues. Of note, KIRC patients with elevated expression of ERVFRD-1 had favorable survival. The multivariate Cox analysis demonstrated that T stage, age, and ERVFRD-1 expression possessed independent prognostic values in KIRC. Collectively, these results demonstrate that ERVFRD-1 holds promise as an attractive and novel prognostic biomarker for KIRC.

To better understand the role of ERVFRD-1 in KIRC, we conducted further investigations into its biological function and signaling pathways. The exploration of differentially expressed genes (DEGs) using GO enrichment and KEGG pathway analyses revealed several immune-related pathways and functions associated with ERVFRD-1. It has been previously demonstrated that epigenetically silenced HERVs can become transcriptionally activated in different cancers, resulting in adaptive immune responses to HERV epitopes (Cherkasova et al., 2013; Downey et al., 2015). Spontaneous T-cell and B-cell responses against HERV antigens are also well-documented (Ebert et al., 2009). Previous studies also suggested ERVFRD-1 as an essential immune regulator, both locally and systemically, through its association with placental exosomes (Lokossou et al., 2020). However, further investigation is required to fully understand the immune regulation of ERVFRD-1 in the context of KIRC.

The tumor microenvironment (TME) is composed of tumor cells and various non-malignant cells, including tumor-infiltrating immune cells (TIICs), fibroblasts, and stromal cells (Wu and Dai, 2017; Zheng et al., 2021). Among these, TIICs have been identified as major regulators of tumor development at different stages (Lei et al., 2020). Building upon our previous findings, we further investigated the relationship between ERVFRD-1 and TIICs in KIRC using GSEA. Our study revealed a positive correlation between ERVFRD-1 expression and mast cells, while a negative correlation was observed with Tregs.

Moreover, previous studies have suggested that the infiltration of dendritic cell resting, mast cell resting, and eosinophilia in KIRC is associated with a favorable prognosis (Pan et al., 2020). Conversely, Tregs are known to be highly immunosuppressive towards effector cells and their increased infiltration during KIRC progression is indicative of a poor prognosis (Şenbabaoğlu et al., 2016; Braun et al., 2021). Additionally, the degree of Treg infiltration is lower in some patients who respond favorably to immune checkpoint inhibitors (ICI) and achieve complete response (Krishna et al., 2021). Taken together, our finding on the prognostic role of ERVFRD-1 supports these previous observations, given its protective effect on overall survival and longer survival outcome.

Checkpoint inhibition therapies targeting PD-L1 or CTLA-4 pathways have shown promising effects in patients with renal cell carcinoma (Choueiri et al., 2016; Motzer et al., 2018). A recent study revealed that responders to these therapies exhibited significantly higher expression of ERV3-2 in tumors compared to non-responders (Panda et al., 2018). Additionally, a TCGA bioinformatics profiling of 18 different cancer types identified specific associations between KIRC and ERV-E, as well as seven different ERV-K families (Rooney et al., 2015). Furthermore, ERV-E env has been identified as a potential target for T-cell-based immunotherapy, specifically in KIRC (Cherkasova et al., 2016). Interestingly, our findings demonstrated a positive correlation between ERVFRD-1 and PD-L1 expression, but not with PD-1 or CTLA-4. Notably, KDR (a gene encoding the VEGF receptor) was found to be the most correlated factor with ERVFRD-1 in our analysis.

In the absence of further supporting evidence, our study aimed to explore the possible relationship between the expression of ERVFRD-1 and immunotherapy in KIRC. Common immunotherapy indices include tumor mutation burden (TMB) and immune checkpoint expression. We aimed to investigate whether ERVFRD-1 has therapeutic value in KIRC immunotherapy by comparing its expression patterns with that of immune checkpoints. Our findings reveal a weak yet statistically significant correlation between TMB and PD-L1 with ERVFRD-1 expression. Next, we sought to determine whether ERVFRD-1 expression has a prognostic value in immunotherapy response in the KIRC population from TCGA dataset. Although our findings did not display a response to PD-L1/PD-1, we must exercise caution in interpreting these results due to our limited sample size, which highlights the need for further validation studies with larger sample sizes. Furthermore, our analysis found that ERVFRD-1 exhibited a higher correlation with other immune stimulators such as CD70 and TNFSF14. These findings suggest that ERVFRD-1 may play a significant role in activating the immune system to achieve an anti-tumor effect in KIRC.

Last but not least, a recent article has reported that the DNA hypomethylating agent decitabine can activate the expression of transposable elements (TEs) such as LINE1, ERV3-2, and ERV4700, as well as antiviral signaling, to enhance the response to immune checkpoint blockade (ICB) in KIRC and primary cells (de Cubas et al., 2020). Thus, targeting ERVs through methylation regulation may represent another potential therapeutic approach for cancer treatment. In light of this, we also investigated the methylation status of ERVFRD-1 in KIRC and found that methylation predominantly drives ERVFRD-1 expression. Additionally, we observed higher levels of methylation in the cg26383454 site were associated with a poor prognosis. Therefore, the use of demethylating agents may also be a viable strategy for improving the prognosis of KIRC patients with high ERVFRD-1 methylation.

5 Conclusions

In summary, this study investigated the expression and potential prognostic value of ERVFRD-1 in patients with KIRC. Through comprehensive analyses of clinical and molecular data from various databases, the study revealed a correlation between ERVFRD-1 expression and immune infiltrates, as well as immune modulators in KIRC. GO and KEGG enrichment analysis revealed that ERVFRD-1 may increase the immune phenotype in KIRC. Furthermore, high ERVFRD-1 expression was associated with a favorable prognostic factor for OS, DSS, and PFI in KIRC patients. Multivariate Cox regression analysis identified ERVFRD-1 expression as an independent prognostic factor for overall survival. A nomogram incorporating ERVFRD-1 expression and other clinical variables was developed to predict patient outcomes, demonstrating moderate predictive accuracy. The findings highlight the potential significance of ERVFRD-1 as a prognostic biomarker and provide a valuable tool for personalized prognosis assessment in KIRC patients.

References

• 1

  AuL.HatipogluE.Robert de MassyM.LitchfieldK.BeattieG.RowanA.et al. (2021). Determinants of anti-PD-1 response and resistance in clear cell renal cell carcinoma. Cancer Cell39 (11), 1497–1518.e1411. doi: 10.1016/j.ccell.2021.10.001

  • CrossRef
  • Google Scholar

• 2

  BaladaE.Ordi-RosJ.Vilardell-TarrésM. (2009). Molecular mechanisms mediated by human endogenous retroviruses (HERVs) in autoimmunity. Rev. Med. Virol.19 (5), 273–286. doi: 10.1002/rmv.622

  • CrossRef
  • Google Scholar

• 3

  BindeaG.MlecnikB.TosoliniM.KirilovskyA.WaldnerM.ObenaufA. C.et al. (2013). Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity39 (4), 782–795. doi: 10.1016/j.immuni.2013.10.003

  • CrossRef
  • Google Scholar

• 4

  BolzeP. A.MommertM.MalletF. (2017). Contribution of syncytins and other endogenous retroviral envelopes to human placenta pathologies. Prog. Mol. Biol. Transl. Sci.145, 111–162. doi: 10.1016/bs.pmbts.2016.12.005

  • CrossRef
  • Google Scholar

• 5

  BraunD. A.StreetK.BurkeK. P.CookmeyerD. L.DenizeT.PedersenC. B.et al. (2021). Progressive immune dysfunction with advancing disease stage in renal cell carcinoma. Cancer Cell39 (5), 632–648.e638. doi: 10.1016/j.ccell.2021.02.013

  • CrossRef
  • Google Scholar

• 6

  CherkasovaE.ScrivaniC.DohS.WeismanQ.TakahashiY.HarashimaN.et al. (2016). Detection of an immunogenic HERV-E envelope with selective expression in clear cell kidney cancer. Cancer Res.76 (8), 2177–2185. doi: 10.1158/0008-5472.can-15-3139

  • CrossRef
  • Google Scholar

• 7

  CherkasovaE.WeismanQ.ChildsR. W. (2013). Endogenous retroviruses as targets for antitumor immunity in renal cell cancer and other tumors. Front. Oncol.3. doi: 10.3389/fonc.2013.00243

  • CrossRef
  • Google Scholar

• 8

  ChinC. H.ChenS. H.WuH. H.HoC. W.KoM. T.LinC. Y. (2014). cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst. Biol.8 Suppl 4 (Suppl 4), S11. doi: 10.1186/1752-0509-8-s4-s11

  • CrossRef
  • Google Scholar

• 9

  ChoueiriT. K.FishmanM. N.EscudierB.McDermottD. F.DrakeC. G.KlugerH.et al. (2016). Immunomodulatory activity of nivolumab in metastatic renal cell carcinoma. Clin. Cancer Res.22 (22), 5461–5471. doi: 10.1158/1078-0432.ccr-15-2839

  • CrossRef
  • Google Scholar

• 10

  ChoueiriT. K.MotzerR. J.RiniB. I.HaanenJ.CampbellM. T.VenugopalB.et al. (2020). Updated efficacy results from the JAVELIN Renal 101 trial: first-line avelumab plus axitinib versus sunitinib in patients with advanced renal cell carcinoma. Ann. Oncol.31 (8), 1030–1039. doi: 10.1016/j.annonc.2020.04.010

  • CrossRef
  • Google Scholar

• 11

  de CubasA. A.DunkerW.ZaninovichA.HongoR. A.BhatiaA.PandaA.et al. (2020). DNA hypomethylation promotes transposable element expression and activation of immune signaling in renal cell cancer. JCI Insight5 (11). doi: 10.1172/jci.insight.137569

  • CrossRef
  • Google Scholar

• 12

  DowneyR. F.SullivanF. J.Wang-JohanningF.AmbsS.GilesF. J.GlynnS. A. (2015). Human endogenous retrovirus K and cancer: Innocent bystander or tumorigenic accomplice? Int. J. Cancer137 (6), 1249–1257. doi: 10.1002/ijc.29003

  • CrossRef
  • Google Scholar

• 13

  EbertP. J.JiangS.XieJ.LiQ. J.DavisM. M. (2009). An endogenous positively selecting peptide enhances mature T cell responses and becomes an autoantigen in the absence of microRNA miR-181a. Nat. Immunol.10 (11), 1162–1169. doi: 10.1038/ni.1797

  • CrossRef
  • Google Scholar

• 14

  GeisF. K.GoffS. P. (2020). Silencing and transcriptional regulation of endogenous retroviruses: an overview. Viruses12 (8). doi: 10.3390/v12080884

  • CrossRef
  • Google Scholar

• 15

  HänzelmannS.CasteloR.GuinneyJ. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf.14, 7. doi: 10.1186/1471-2105-14-7

  • CrossRef
  • Google Scholar

• 16

  JiangP.GuS.PanD.FuJ.SahuA.HuX.et al. (2018). Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat. Med.24 (10), 1550–1558. doi: 10.1038/s41591-018-0136-1

  • CrossRef
  • Google Scholar

• 17

  KimM. C.JinZ.KolbR.BorcherdingN.ChatzkelJ. A.FalzaranoS. M.et al. (2021). Updates on immunotherapy and immune landscape in renal clear cell carcinoma. Cancers (Basel)13 (22). doi: 10.3390/cancers13225856

  • CrossRef
  • Google Scholar

• 18

  KotechaR. R.MotzerR. J.VossM. H. (2019). Towards individualized therapy for metastatic renal cell carcinoma. Nat. Rev. Clin. Oncol.16 (10), 621–633. doi: 10.1038/s41571-019-0209-1

  • CrossRef
  • Google Scholar

• 19

  KrishnaC.DiNataleR. G.KuoF.SrivastavaR. M.VuongL.ChowellD.et al. (2021). Single-cell sequencing links multiregional immune landscapes and tissue-resident T cells in ccRCC to tumor topology and therapy efficacy. Cancer Cell39 (5), 662–677.e666. doi: 10.1016/j.ccell.2021.03.007

  • CrossRef
  • Google Scholar

• 20

  LeiX.LeiY.LiJ. K.DuW. X.LiR. G.YangJ.et al. (2020). Immune cells within the tumor microenvironment: Biological functions and roles in cancer immunotherapy. Cancer Lett.470, 126–133. doi: 10.1016/j.canlet.2019.11.009

  • CrossRef
  • Google Scholar

• 21

  LiJ. Y.ChenY. P.LiY. Q.LiuN.MaJ. (2021). Chemotherapeutic and targeted agents can modulate the tumor microenvironment and increase the efficacy of immune checkpoint blockades. Mol. Cancer20 (1), 27. doi: 10.1186/s12943-021-01317-7

  • CrossRef
  • Google Scholar

• 22

  LokossouA. G.ToudicC.NguyenP. T.ElisseeffX.VargasA.RassartÉ.et al. (2020). Endogenous retrovirus-encoded Syncytin-2 contributes to exosome-mediated immunosuppression of T cells†. Biol. Reprod.102 (1), 185–198. doi: 10.1093/biolre/ioz124

  • CrossRef
  • Google Scholar

• 23

  LotiaS.MontojoJ.DongY.BaderG. D.PicoA. R. (2013). Cytoscape app store. Bioinformatics29 (10), 1350–1351. doi: 10.1093/bioinformatics/btt138

  • CrossRef
  • Google Scholar

• 24

  LoveM. I.HuberW.AndersS. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol.15 (12), 550. doi: 10.1186/s13059-014-0550-8

  • CrossRef
  • Google Scholar

• 25

  MakhovP.JoshiS.GhataliaP.KutikovA.UzzoR. G.KolenkoV. M. (2018). Resistance to systemic therapies in clear cell renal cell carcinoma: mechanisms and management strategies. Mol. Cancer Ther.17 (7), 1355–1364. doi: 10.1158/1535-7163.mct-17-1299

  • CrossRef
  • Google Scholar

• 26

  MangeneyM.RenardM.Schlecht-LoufG.BouallagaI.HeidmannO.LetzelterC.et al. (2007). Placental syncytins: Genetic disjunction between the fusogenic and immunosuppressive activity of retroviral envelope proteins. Proc. Natl. Acad. Sci. U.S.A.104 (51), 20534–20539. doi: 10.1073/pnas.0707873105

  • CrossRef
  • Google Scholar

• 27

  MochH.AminM. B.BerneyD. M.CompératE. M.GillA. J.HartmannA.et al. (2022). The 2022 world health organization classification of tumours of the urinary system and male genital organs-part A: renal, penile, and testicular tumours. Eur. Urol82 (5), 458–468. doi: 10.1016/j.eururo.2022.06.016

  • CrossRef
  • Google Scholar

• 28

  ModhukurV.IljasenkoT.MetsaluT.LokkK.Laisk-PodarT.ViloJ. (2018). MethSurv: a web tool to perform multivariable survival analysis using DNA methylation data. Epigenomics10 (3), 277–288. doi: 10.2217/epi-2017-0118

  • CrossRef
  • Google Scholar

• 29

  MotzerR. J.TannirN. M.McDermottD. F.Arén FronteraO.MelicharB.ChoueiriT. K.et al. (2018). Nivolumab plus Ipilimumab versus Sunitinib in Advanced Renal-Cell Carcinoma. N Engl. J. Med.378 (14), 1277–1290. doi: 10.1056/NEJMoa1712126

  • CrossRef
  • Google Scholar

• 30

  NairV. P.LiuH.CiceriG.JungverdorbenJ.FrishmanG.TchieuJ.et al. (2021). Activation of HERV-K(HML-2) disrupts cortical patterning and neuronal differentiation by increasing NTRK3. Cell Stem Cell28 (9), 1671–1673. doi: 10.1016/j.stem.2021.05.003

  • CrossRef
  • Google Scholar

• 31

  NambaS.UenoT.KojimaS.KobayashiK.KawaseK.TanakaY.et al. (2021). Transcript-targeted analysis reveals isoform alterations and double-hop fusions in breast cancer. Commun. Biol.4 (1), 1320. doi: 10.1038/s42003-021-02833-4

  • CrossRef
  • Google Scholar

• 32

  NgK. W.BoumelhaJ.EnfieldK. S. S.AlmagroJ.ChaH.PichO.et al. (2023). Antibodies against endogenous retroviruses promote lung cancer immunotherapy. Nature616 (7957), 563–573. doi: 10.1038/s41586-023-05771-9

  • CrossRef
  • Google Scholar

• 33

  OkaharaG.MatsubaraS.OdaT.SugimotoJ.JinnoY.KanayaF. (2004). Expression analyses of human endogenous retroviruses (HERVs): tissue-specific and developmental stage-dependent expression of HERVs. Genomics84 (6), 982–990. doi: 10.1016/j.ygeno.2004.09.004

  • CrossRef
  • Google Scholar

• 34

  PanQ.WangL.ChaiS.ZhangH.LiB. (2020). The immune infiltration in clear cell Renal Cell Carcinoma and their clinical implications: A study based on TCGA and GEO databases. J. Cancer11 (11), 3207–3215. doi: 10.7150/jca.37285

  • CrossRef
  • Google Scholar

• 35

  PandaA.de CubasA. A.SteinM.RiedlingerG.KraJ.MayerT.et al. (2018). Endogenous retrovirus expression is associated with response to immune checkpoint blockade in clear cell renal cell carcinoma. JCI Insight3 (16). doi: 10.1172/jci.insight.121522

  • CrossRef
  • Google Scholar

• 36

  PortaC.CosmaiL.LeibovichB. C.PowlesT.GallieniM.BexA. (2019). The adjuvant treatment of kidney cancer: a multidisciplinary outlook. Nat. Rev. Nephrol.15 (7), 423–433. doi: 10.1038/s41581-019-0131-x

  • CrossRef
  • Google Scholar

• 37

  PowlesT.TomczakP.ParkS. H.VenugopalB.FergusonT.SymeonidesS. N.et al. (2022). Pembrolizumab versus placebo as post-nephrectomy adjuvant therapy for clear cell renal cell carcinoma (KEYNOTE-564): 30-month follow-up analysis of a multicentre, randomised, double-blind, placebo-controlled, phase 3 trial. Lancet Oncol.23 (9), 1133–1144. doi: 10.1016/s1470-2045(22)00487-9

  • CrossRef
  • Google Scholar

• 38

  RangelS. C.da SilvaM. D.da SilvaA. L.Dos SantosJ. M. B.NevesL. M.PedrosaA.et al. (2022). Human endogenous retroviruses and the inflammatory response: A vicious circle associated with health and illness. Front. Immunol.13. doi: 10.3389/fimmu.2022.1057791

  • CrossRef
  • Google Scholar

• 39

  RobinsonM. D.McCarthyD. J.SmythG. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics26 (1), 139–140. doi: 10.1093/bioinformatics/btp616

  • CrossRef
  • Google Scholar

• 40

  RooneyM. S.ShuklaS. A.WuC. J.GetzG.HacohenN. (2015). Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell160 (1-2), 48–61. doi: 10.1016/j.cell.2014.12.033

  • CrossRef
  • Google Scholar

• 41

  ŞenbabaoğluY.GejmanR. S.WinerA. G.LiuM.Van AllenE. M.de VelascoG.et al. (2016). Tumor immune microenvironment characterization in clear cell renal cell carcinoma identifies prognostic and immunotherapeutically relevant messenger RNA signatures. Genome Biol.17 (1), 231. doi: 10.1186/s13059-016-1092-z

  • CrossRef
  • Google Scholar

• 42

  SubramanianA.TamayoP.MoothaV. K.MukherjeeS.EbertB. L.GilletteM. A.et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U.S.A.102 (43), 15545–15550. doi: 10.1073/pnas.0506580102

  • CrossRef
  • Google Scholar

• 43

  SungH.FerlayJ.SiegelR. L.LaversanneM.SoerjomataramI.JemalA.et al. (2021). Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J. Clin.71 (3), 209–249. doi: 10.3322/caac.21660

  • CrossRef
  • Google Scholar

• 44

  SuntsovaM.GarazhaA.IvanovaA.KaminskyD.ZhavoronkovA.BuzdinA. (2015). Molecular functions of human endogenous retroviruses in health and disease. Cell Mol. Life Sci.72 (19), 3653–3675. doi: 10.1007/s00018-015-1947-6

  • CrossRef
  • Google Scholar

• 45

  TurnerG.BarbulescuM.SuM.Jensen-SeamanM. I.KiddK. K.LenzJ. (2001). Insertional polymorphisms of full-length endogenous retroviruses in humans. Curr. Biol.11 (19), 1531–1535. doi: 10.1016/s0960-9822(01)00455-9

  • CrossRef
  • Google Scholar

• 46

  WeyererV.StrisselP. L.StöhrC.EcksteinM.WachS.TaubertH.et al. (2021). Endogenous retroviral-K envelope is a novel tumor antigen and prognostic indicator of renal cell carcinoma. Front. Oncol.11. doi: 10.3389/fonc.2021.657187

  • CrossRef
  • Google Scholar

• 47

  WuT.DaiY. (2017). Tumor microenvironment and therapeutic response. Cancer Lett.387, 61–68. doi: 10.1016/j.canlet.2016.01.043

  • CrossRef
  • Google Scholar

• 48

  YoshiharaK.ShahmoradgoliM.MartínezE.VegesnaR.KimH.Torres-GarciaW.et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun.4, 2612. doi: 10.1038/ncomms3612

  • CrossRef
  • Google Scholar

• 49

  YuH.LiuT.ZhaoZ.ChenY.ZengJ.LiuS.et al. (2014). Mutations in 3'-long terminal repeat of HERV-W family in chromosome 7 upregulate syncytin-1 expression in urothelial cell carcinoma of the bladder through interacting with c-Myb. Oncogene33 (30), 3947–3958. doi: 10.1038/onc.2013.366

  • CrossRef
  • Google Scholar

• 50

  YuG.WangL. G.HanY.HeQ. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. Omics16 (5), 284–287. doi: 10.1089/omi.2011.0118

  • CrossRef
  • Google Scholar

• 51

  ZhangM.LiangJ. Q.ZhengS. (2019). Expressional activation and functional roles of human endogenous retroviruses in cancers. Rev. Med. Virol.29 (2), e2025. doi: 10.1002/rmv.2025

  • CrossRef
  • Google Scholar

• 52

  ZhengB.XieF.ChengF.WangJ.YaoZ.HeW.et al. (2021). Integrative analysis of immune-related genes in the tumor microenvironment of renal clear cell carcinoma and renal papillary cell carcinoma. Front. Mol. Biosci.8. doi: 10.3389/fmolb.2021.760031

  • CrossRef
  • Google Scholar