Locus-Specific Characterization of Human Endogenous Retrovirus Expression in Prostate, Breast, and Colon Cancers

Prostate, breast, and colon cancers together account for over 500,000 new cancer cases in the United States each year (1). Recently, human endogenous retroviruses (HERV), a class of transposable elements resulting from the integration of ancient exogenous retroviruses into the human genome (2), have been implicated in cancer mechanisms. HERVs comprise approximately 8% to 9% of the human genome and have been incorporated into primate genomes via retrotransposition for at least 40 million years (2). Sequence degradation over time has rendered the majority of protein-coding HERV sequences inactive; however, HERV transcription has been identified in all commonly studied human tissue types (3). While HERV nomenclature varies in the literature, HERVs separate phylogenetically into over 40 distinct groups, the most commonly studied being the HERV-K group and its subgroup HML-2.

Previous in vitro work in tumor cell lines and tissue samples has identified evidence of HERV transcription in prostate, breast, and colon cancers. Individual HERV-K elements including HERV-K102 and HERV-K118, as well as solo long terminal repeats (LTR), are transcribed in prostate cancer cell lines (4). Expression of an antigen related to the HERV-K gag gene has also been found in prostate tumors (5). In addition, HERVK expression has been found in other urologic cancers (6, 7), and the HERV-K(HML-2) protein Rec has been suggested to be oncogenic in these cancers by derepressing the androgen receptor (8). In breast cancers, HERV transcription has been found in nonmalignant breast tissue, as well as lower overall HERV transcription in tumor tissues, with the latter being heterogeneous across patients (9). Further studies have identified evidence of HERV-K transcripts including HERV-K102 and HERV-K(HML-2) loci 1q22 and 3q21.2 (10–12) in breast cancers, as well as evidence of HERV-K env transcription (13, 14). In colorectal cancers, several HERV-H elements are expressed in both tumor and normal tissues, particularly at the locus Xp22.32 (15, 16), as well as the ERV3–1 Env protein (17).

Unveiling HERV expression patterns across multiple cancer types is important to understanding the broader role of these elements in cancer progression. Toward this aim, large biomedical databases have been utilized to conduct pan-cancer studies of transposable element expression (18–21). However, identifying locus-specific HERV expression in transcriptomics data remains challenging from a bioinformatics standpoint due to high sequence similarity among HERV groups, and therefore, the vast majority of these studies report group-level expression patterns. As such, novel computational frameworks are needed in tandem with large -omics datasets to identify locus-specific HERV expression patterns across human cancers.

While HERV transcription signatures are useful as prognostic biomarkers (22–24), elucidating the underlying mechanisms of HERV activation and function in the context of cancer development is also essential to treatment applications. As methylation is a key mechanism in silencing HERV function, hypomethylation patterns observed in tumor cells may induce HERV expression. HERV expression has also been shown to contribute to immunomodulation of host cells and HERVs have been implicated in the regulation of innate immunity (25, 26). Furthermore, a viral mimicry phenotype has been identified in subsets of metastatic tumors, which correspond to increased HERV expression, further suggesting the contribution of HERVs to inflammatory immune responses supporting tumor survival (27). Thus, an understanding of host cell factors including epigenetics and immunogenicity, in addition to HERV expression, is highly relevant to translational applications.

Here, we utilize Telescope (28) to present a locus-specific characterization of differential HERV expression in paired tumor-normal prostate, breast, and colon cancer RNA-sequencing (RNA-seq) data from The Cancer Genome Atlas (TCGA; ref. 29). We then combine these results with differential expression analysis of host genes to identify potential host genomic targets of HERV expression and regulatory mechanisms. Finally, we analyze broader patterns in dysregulated HERV elements across these three cancer types to identify potential mechanisms of HERV activation associated with cancers.

Paired-end RNA-seq data (Illumina) were downloaded in BAM format (*.bam) from TCGA (29) using the Genomic Data Commons (GDC) Data Transfer Tool (30). In total, RNA-seq files from paired primary tumor (non-metastatic) and solid tissue normal (noncancerous tissues adjacent to the tumor) samples from 52 prostate cancer (TCGA-PRAD) cases, 111 breast cancer (TCGA-BRCA) cases, and 24 colon cancer (TCGA-COAD) cases were pulled. FFPE and metastatic samples were excluded from further analysis, and in cases where multiple tumor files remained for a single case, per-base sequence quality (evaluated using FastQC v0.11.9; ref. 31) was used to select one file for use in this study. Corresponding clinical and demographic metadata were downloaded using TCGABiolinks v3.11 (32) and all related analyses were performed in R v4.0.0 (33) (Supplementary Table S1).

Downloaded BAM files were converted to FASTQ format using SAMtools v. 1.10 (34). FASTQ files were then trimmed using Trimmomatic v. 0.39 (35) with a leading/trailing threshold of 3, sliding window of 4bp with a threshold of 14, minimum read length of 36 bp, and Illumina adapters removed. In TCGA-PRAD data, the percentage of reads remaining after trimming ranged from 71.30% to 97.50% (average 95.15 % ± 2.71) and the number of trimmed reads per file ranged from 24,556,156 to 112,956,172 (average 72,099,135 ± 12,007,080). In TCGA-BRCA data, the percentage of reads remaining after trimming ranged from 65.02 to 97.08% (average 93.41% ± 3.24) and the number of trimmed reads per file ranged from 28,327,165 to 150,387,689 (average 70,310,030 ± 17,254,417). In TCGA-COAD data, the percentage of reads remaining after trimming ranged from 88.44 to 96.85% (average 94.02% ± 2.57) and the number of trimmed reads per file ranged from 28,901,359 to 101,130,017 (average 56,122,694.8 ± 13,567,413.5).

Trimmed reads were aligned to Hg38 using Bowtie2 v2.3.5.1 (36) with a very sensitive local alignment search, up to 100 alignments per read, and with fragments having at least 95% identity reported. Alignment rates averaged 92.9% (±0.94%) for TCGA-PRAD, 96.0% (±0.025%) for TCGA-BRCA, and 92.1% (±0.72%) for TCGA-COAD. The Telescope v1.0.3–0 (28) assign module was then used for transposable element quantification using the HERV annotation described in ref. 28 and a theta prior of 200,000. The Telescope algorithm is unique in identifying locus-specific HERV expression, as opposed to groups of HERVs, as is common in other studies. In addition, Telescope uses a novel Bayesian reassignment method to map fragments onto transcripts, allowing for the utilization of more data than other approaches, which either are limited to mapping just unique reads, thereby, using a very limited dataset or incorporate all reads, but without the reassignment step, resulting in low resolution to differentiate among alternative mapping options.

Kallisto v0.46.2 (37) was used to quantify gene expression data using default parameters and the Ensembl v96 human transcriptome assembly indexed with default settings. Tximport v1.16.1 (38) was used to read abundance files into R. Average pseudoalignment rate was 92.2% (±2.3%) for TCGA-PRAD data, 90.5% (±5.3%) for TCGA-BRCA data, and 90.9% (±7.5%) for TCGA-COAD data.

Differential expression analysis of retrotranscriptomic and transcriptomic data between tumor and normal samples was performed using DESeq2 v1.28.1 (39) using a negative binomial model and Wald test with a significance cutoff of P = 0.05, using the patient identifier and condition (tumor or normal tissue) as factors in the DESeq model (∼ patient + condition). Variance stabilizing transformation was used prior to principal component analysis (PCA) and hierarchical clustering. The R packages ggplot2 v3.2.2 (40) and EnhancedVolcano v1.6.0 (41) were used in addition to DESeq2 functions to visualize PCA and volcano plots, respectively, and the circlize v0.4.10 (42) and GenomicRanges v1.40.0 (43) packages were used to construct circos plots of differential gene and HERV expression patterns. Lastly, the pheatmap v1.0.12 (44) package was used to visualize sample-wise expression and perform hierarchical clustering for top differentially expressed HERVs.

Gene set enrichment analyses (GSEA) of differentially expressed HERVs and differentially expressed genes were performed using GSEAPreranked as implemented in GSEA v4.0.3 (45). HERVs and genes, as applicable, were ranked by sign(log₂(fold change))*−log₁₀(P_adj), using fold change metrics and P values from DESeq2. For HERVs, a custom gene set was generated from the HERV annotation described in ref. 28. For genes, gene ontology sets from the Molecular Signatures Database (MSigDB; refs. 45, 46), as provided by GSEA, were used.

Survival analysis was performed using the survival R package v. 3.2.3 (47) using the log-rank test for by-group comparisons. Kaplan–Meier survival curves were visualized using ggfortify v. 0.4.10 (48). Survival times were censored at 5 years, with patients listed as alive censored at the date of the last follow-up.

Annotations of the intersected and nearby genes for the transposable element reference utilized here were constructed using the ENSEMBL Hg38 release 99. Annotations of enhancer regions were constructed using FANTOM5 (49).

In total, out of 4,965 HERVs identified, 720 were differentially expressed between tumor and normal prostate cancer samples, with 361 upregulated in tumor tissue and 259 downregulated (Fig. 1A; Table 1). A full list of differentially expressed HERVs is available in the Supplementary Data. The HERV groups MER4, HERVL/ERVL, HERVK, and HERVH were most prevalent among both up- and downregulated HERV elements (Supplementary Data; Supplementary Table S2), however, only HERVE was significantly enriched in downregulated HERVs (q < 0.25), and no groups were significantly enriched in upregulated HERVs (Supplementary Data; Supplementary Table S3). PCA of the 500 most variable HERVs shows separation of tumor and normal HERV expression data in the majority of samples, and also a group of seven normal tissue samples that cluster separately (Fig. 1B). These seven samples additionally cluster together in hierarchical clustering analysis of the top 100 most significantly differentially expressed HERV loci (Supplementary Data; Supplementary Fig. S1), which shows notable upregulation of HML3 10q11.21a, HARLEQUIN 1q32.1, HERVIP10F 19p13.2b, and PABLA 14q11.2 and downregulation of HERVL18 2q31.1 in these samples compared with other normal tissue samples. Secondary Gleason grade, age at diagnosis, and race were not significantly different between cases corresponding to each group [χ² test, P values in (0.21 to 0.79)]. Survival data were not available for this dataset.

In addition, 12,861 genes were identified as differentially expressed in tumor samples as compared with normal samples, with 6,423 being upregulated in tumor samples and 6,438 downregulated (Fig. 1C; Supplementary Data; Supplementary Fig. S2A). Gene ontology enrichment analysis identified the enrichment of multiple processes relevant to HERV expression in upregulated genes, including viral protein expression (e.g., GO:0019080, q < 0.05), epigenetic gene regulation (e.g., GO:0016458, GO:0045814, GO:0040029; q < 0.05), and immune response, including innate immune response (e.g., GO:0002251, GO:0002227; q < 0.25; Supplementary Data). Genes related to immunoglobulin secretion were also enriched in downregulated genes (e.g., q < 0.25), as well as interferon response (e.g., GO:0035455, q < 0.25) and inflammation (GO:0002269, q < 0.05). In addition, two TET genes (TET1 and TET3), which are associated with demethylation and are frequently mutated in cancer cells, were significantly upregulated in tumor samples. Several APOBEC genes, which are also involved in DNA demethylation pathways, were also differentially expressed, including upregulation of APOBEC3B and downregulation of APOBEC3C, APOBEC3G, APOBEC3F, APOBEC3D, APOBEC3H, and APOBEC2 in tumor samples (P < 0.05). Several of these genes are also involved in RNA editing and can act as viral restriction factors. Furthermore, the antiviral signaling gene MAVS, interferon regulating factor IFN7, and two interferon genes (IFNAR1 and IFNGR2) were upregulated (P < 0.05). The seven normal tissue samples that clustered separately by HERV expression also clustered separately by host gene expression in PCA and hierarchical clustering (Supplementary Data; Supplementary Fig. S2B and S3). AOX1, which has been associated with prostate cancer survival times, displays higher expression levels in this subgroup, as well as GSTP1, the methylation of which has been associated with prostate tumor carcinogenesis. GCNT4, a colon cancer biomarker, additionally shows higher expression. Furthermore, the two most significantly upregulated HERVs in PRAD—both of which are HERV-K elements—are located on the q arm of chromosome 22, which is also the site of a TMPRSS2 and ERG gene fusion leading to overexpression of ERG, an oncogene. This fusion is among the most common in prostate cancer, and ERG was additionally found to be significantly overexpressed in tumor samples (P < 0.001).

We then used HERV meta-annotations provided by Telescope to examine functional properties of differentially expressed elements, as well as nearby or intersected genes. Of 720 differentially expressed HERVs, 262 are located in exonic regions, 232 are located in intronic regions, and 226 are located in intergenic regions. A total of 67 are located in enhancer regions. Furthermore, these HERVs intersect 271 protein-coding genes and 242 long noncoding RNAs (lncRNA). Several of these regions correspond to known gene expression patterns in prostate and other cancers. Specifically, among upregulated HERVs, HERVK11 22q11.21 intersects TXNRD2, a protein-coding gene, which was also significantly upregulated in tumor samples (P < 0.001). TXNRD2 expression has been linked with tumor growth and suggested as a drug target. In addition, HML2 22a11.23 intersects PCAT14, a lncRNA that has been identified as a prostate cancer biomarker. HML1 16p12.3a intersects the genes ACSM3 and ACSM1, both of which have been suggested as clinical biomarkers for prostate cancer. ACSM1 was significantly upregulated in tumor samples (P < 0.001). Among downregulated HERVs, PABLA 14q11.2 intersects AL161668.3, a lncRNA, which has been shown to be downregulated in head and neck cancer. HARLEQUIN 10q23.1 intersects CERNA2, a lncRNA associated with late-stage and metastatic cervical cancers and malignant gastric cancers. In total, 155 differentially expressed HERVs intersected a gene that was also differentially expressed, with 142 of 155 being differentially expressed in the same direction. Furthermore, several differentially expressed HERVs are proximal to relevant genes. Notably, PABLB 4p15.2, which was upregulated in tumor samples (P < 0.001), is 15.3 kb 3′ of DHX15, the expression of which is associated with androgen receptor activation and prostate cancer progression. In addition, HML3 8p11.22, also upregulated (P < 0.001) is 26.3 kb 5′ of IDO1, the upregulation of which is implicated in a mechanism for epithelial-to-mesenchymal transition in prostate cancer cells. In total, there were 75 differentially expressed HERVs for which the closest 5′ gene was also differentially expressed, and there were 77 differentially expressed HERVs for which the closest 3′ gene was also differentially expressed.

In total, of 6,330 HERVs identified, 1,547 HERVs were differentially expressed in tumor versus normal TCGA-BRCA samples (P < 0.05), including 640 upregulated and 907 downregulated HERVs (Table 2; Fig. 2A). A full list of differentially expressed HERVs is available in the supplemental data. MER4, HERVL/ERVL, HERVH, and HERVK were most prevalent among both up- and downregulated HERVs (Supplementary Data; Supplementary Table S4). GSEA results indicate that groups HERVL and ERV1 are enriched among upregulated HERVs, while HERVH is enriched in downregulated HERVs (false discovery rate q < 0.25; Supplementary Data; Supplementary Table S5). PCA of the 500 most variable HERVs shows separation of tumor and normal samples, with some separation of normal samples along principal component 2 (Fig. 2B). This separation corresponds with the clustering of 14 normal samples separated from others in hierarchical clustering of the 100 most variable HERVs (Supplementary Data; Supplementary Fig. S4), in which decreased expression of several HERVs relative to other normal samples is visible. Survival analysis revealed no significant difference in survival times of these 14 cases versus others (P = 0.1; Supplementary Data; Supplementary Fig. S5). Tumor stage, menopause stage, age at diagnosis, and race were not significantly different between cases corresponding to each group [χ² test, P values in (0.17–1)]. However, a greater proportion of patients in the smaller subset had received hormone therapy as a part of treatment (χ² test, P = 0.03).

In total, 18,425 genes were identified as differentially expressed in the TCGA-BRCA data, including 8,609 upregulated and 9,816 downregulated in primary tumors (Fig. 2C; Supplementary Data; Supplementary Fig. S6A). Similar to the prostate cancer data, pathways related to epigenetic gene regulation (e.g., GO:0016458, GO:0045815, GO:0045814, GO:0040029; q < 0.05) and innate immune response (e.g., GO:0045089, GO:0002218; q < 0.25) were significantly enriched in upregulated genes, as is consistent with proposed mechanisms of HERV activation (Supplementary Data). BRCA1 and BRCA2 were both significantly upregulated in tumor samples. Several methylation-associated genes were also dysregulated, including upregulation of TET3, APOBEC3B, APOBEC3A, APOBEC3H, and APOBEC4 as well as downregulation of APOBEC3C, APOBEC2, and APOBEC3G (P < 0.05). In addition, IRF7 and several interferon genes (IFN1, IFNL3P1, IFNGR2, and IFNG) were overexpressed (P < 0.05). The clustering of normal samples in PCA is consistent with that found in HERV expression PCA, and a majority of the 14 samples in the separate subgroup cluster together in hierarchical clustering (Supplementary Data; Supplementary Fig. S6B and S7). The samples in the subgroup and those that cluster with them display a notable increase in LEP expression, which is known to be procarcinogenic in breast cancers.

Of the 1,547 differentially expressed HERVs identified, 497, 554, and 496 are located in exonic, intronic, and intergenic regions, respectively. In addition, 133 are located in enhancer regions. The identified HERVs intersect 651 protein coding genes and 443 lncRNAs. Of note in upregulated HERVs, HARLEQUIN 17q21.31 intersects BRCA1. LTR25 7q21.12 and MER61 7q21.12 both intersect TP53TG1, a lncRNA and component of the p53 pathway, which inhibits tumor growth and is inactivated in human cancer cells via methylation. PABLA 2p35.1 intersects RRM2, which promotes breast cancer metastasis and was upregulated in tumor samples (P < 0.001). In addition, HERV3 21q22.13 intersects KCNJ15, which has been associated with esophageal and renal cancers. KCNJ15 was also upregulated in tumor samples (P < 0.001). ERV316A3 8q23.1c, which was downregulated in tumor samples, intersects PKHD1L1, which is commonly mutated in triple-negative breast cancers. PKHD1L1 was also downregulated in tumor samples (P < 0.001). In total, there were 387 differentially expressed HERVs for which the intersected gene was also differentially expressed, with 336 of 387 being differentially expressed in the same direction. HERVH 6p25.1, which was downregulated in tumor samples, is 20.2 kb 5′ of NRN1, which is known to be hypermethylated in breast cancer and was also downregulated in tumor samples (P < 0.001). In total, there were 84 differentially expressed HERVs for which the closest 5′ gene was also differentially expressed and 73 differentially expressed HERVs for which the closest 3′ gene was also differentially expressed.

In the TCGA-COAD data, of 3,586 HERVs identified, 575 were differentially expressed in tumor versus normal samples, with 283 upregulated and 292 downregulated (P < 0.05; Table 3; Fig. 3A). A full list of differentially expressed HERVs is available in the Supplementary Data. The HERV groups MER4, HERVL/ERVL, HERVH, and HERVK were most prevalent among both up- and downregulated HERVs (Supplementary Data; Supplementary Table S6). GSEA results indicate enrichment of PAB and ERV1 groups in upregulated HERVs (Supplementary Data; Supplementary Table S7). PCA of HERV expression data shows a clear distinction between tumor and normal samples, with additional separation of two tumor samples along PC2 (Fig. 3B). These two samples, along with two others, cluster with normal samples in hierarchical clustering, although separation between tumor and normal samples is otherwise distinct (Supplementary Data; Supplementary Fig. S8). Overall differences in HERV expression data are less apparent than in PRAD or BRCA datasets, although notable overexpression of two HERVH elements (Xp22.32a and 20p11.23b) is seen in several tumor samples.

In addition, 12,701 genes were differentially expressed in the TCGA-COAD data, including 5,851 upregulated and 6,850 downregulated genes (Fig. 3C; Supplementary Data; Supplementary Fig. S9A). As before, upregulated genes are enriched for pathways relevant to epigenetic gene regulation (e.g., GO:0016458, GO:0045815, GO:0045814, GO:0040029; q < 0.25). Notably, four of the top ten most significantly enriched pathways in downregulated genes are related to immune response and immunoglobulin circulation (GO:0019814, GO:0034987, GO:0042571, GO:0002455; q < 0.05; Supplementary Data). TET2 was significantly downregulated in tumor samples (P < 0.001), which would promote DNA methylation, as well as APOBEC3C, APOBEC2, and APOBEC3G (P < 0.05). Several interferon genes (IFNWP19, IFNLP31, and IFNG) and interferon regulatory factor IRF7 were also significantly overexpressed (P < 0.05). Hierarchical clustering and PCA of gene expression data show clear separation of tumor and normal samples, as well as distinct differences in gene expression patterns (Supplementary Data; Supplementary Fig. S9B and S10).

The 575 differentially expressed HERVs are distributed across exonic, intronic, and intergenic regions with 218, 172, and 185 per region, respectively. Of these, 63 are located in enhancer regions. These HERVs intersect 231 protein-coding genes and 178 lncRNAs. Included among these, as in the breast cancer data, is MER61 7q21.12, which intersects TP53TG1. In addition, HERVH 4q22.1c, which was downregulated in tumor samples, intersects ABCG2, which was also significantly downregulated (P < 0.001) and has been linked to treatment resistance in colon cancer (50, 51). HERVH 4q22.3a, which was also downregulated in tumor samples, intersects the downregulated gene (P < 0.001) HPGDS, which has been shown to inhibit intestinal tumors in a mouse model. In total, there were 119 differentially expressed HERVs for which the intersected gene was also differentially expressed, with 114 of 119 being differentially expressed in the same direction. In addition, there were 70 differentially expressed HERVs for which the closest 5′ gene was also differentially expressed and 58 for which the closest 3′ gene was also differentially expressed. Included among these is HERVL 4q32.1a, which was downregulated in tumor samples and is 1.1 kb 5′ of TDO2, which was upregulated in tumor samples (P < 0.01) and has been associated with advanced stage colorectal cancers.

Many HERV loci were identified as differentially expressed in more than one of prostate, breast, and colon cancers, including 155 loci that were differentially expressed in all three cancer types (Fig. 4A). Of these, 114 were differentially expressed in the same direction in all three cancer types. The majority were of the groups MER4 (26/155), HERVH (25/155), ERVL (18/155), HERVK (16/155), HERVL (13/155), and ERV3 (11/155; Fig. 4B). Of these HERV loci, 72 are located in exonic regions, 43 are located in intronic regions, and 40 are located in intergenic regions. These elements intersect 77 protein-coding genes and 53 lncRNAs. Notably, one HERV that is differentially expressed in all three cancer types (HARLEQUIN 1q32.1) intersects the tumor suppressor gene BRCA1. This element was upregulated in PRAD and COAD samples, but downregulated in BRCA samples. Furthermore, 279 HERV loci were differentially expressed in both PRAD and BRCA datasets, 183 were differentially expressed in both BRCA and COAD datasets, and 58 were differentially expressed in both COAD and PRAD datasets. Patterns in the prevalence of HERV groups among up- and downregulated HERVs are similar across all three cancer types, with MER4, HERVL, HERVK, HERVH, and ERVL being the five most frequently dysregulated groups in all datasets.

Here, we present locus-specific differential HERV expression patterns in prostate, breast, and colon cancers. In general, MER4, HERVK, and HERVL/ERVL elements were most commonly dysregulated in tumor tissues across breast, prostate, and colon cancers, with over 600 unique elements being dysregulated in more than one cancer type. In addition, we identified host gene expression and analyzed the intersection sites of dysregulated HERV loci in order to address potential pathogenic mechanisms.

In the TCGA-PRAD dataset, a total of 720 differentially expressed HERVs were identified, as well as 12,861 host genes. A majority of the differentially expressed HERVs were upregulated versus downregulated elements (361 and 359, respectively), as has been shown to be common in many cancer types (19). Gene-set enrichment analysis found that no HERV groups were enriched in upregulated HERVs and only HERVE was enriched in downregulated HERVs (q < 0.25), likely owing to the prevalence of more frequently dysregulated HERV groups (e.g., MER4, HERVL/ERVL, HERVK, and HERVH) in both up- and downregulated loci. Two distinct clusters of normal samples were observed in both HERV and host gene expression data, which did not correlate with sequencing batches or clinical variables. This clustering was consistent in PCA and hierarchical clustering of both HERVs and host genes. In addition, several cancer-associated genes show overexpression in the smaller (n = 7) subset, including AOX1, GSTP1, and GCNT4, with these findings together suggesting that the clustering pattern has biological significance. This observation is consistent with previous empirical studies that found subgroups of patients with significantly different HERV expression profiles, including in normal tissue samples (9), and those that have identified prognostic markers in tumor-adjacent normal tissues (52). Furthermore, another recent study of TCGA data identified heterogeneity in HERV expression across patients, as is consistent with this observation (21). The correspondence between HERV and gene expression profiles (Fig. 1C)—as well as the proximity to and/or intersection of differentially expressed, cancer-associated genes by the identified HERVs (152 and 155 unique HERVs, respectively)—further suggests a functional link between HERV expression and cancer pathways. A majority (513/720) of differentially expressed HERVs intersect either protein-coding genes or lncRNAs; the latter being typically associated with HERVs and activated in cancers (53). Several known prostate cancer-associated genes and lncRNAs were intersected by or located immediately 5′ or 3′ to a differentially expressed HERV, including TXNRD2, PCAT14, ACSM1/3, AL161668.3, CERNA2, DHX15, and IDO1.

In the TCGA-BRCA dataset, a total of 1,547 differentially expressed HERVs were identified, as well as 18,425 genes. The most significantly differentially expressed HERVs are predominantly downregulated (907/1,546; Fig. 2B), as is consistent with previous studies (22, 31). HERVL and ERV1 were enriched in upregulated HERVs, while HERVH was enriched in downregulated HERVs (q < 0.25). As in the TCGA-PRAD data, a subgroup (n = 14) of normal samples was identified by PCA and hierarchical clustering that was consistent between HERV and gene expression data. The presence of a subset of samples with different HERV expression profiles is consistent with a previous empirical study in breast cancer cell lines (9). Here, this subgroup also displays a clear increase in LEP expression, which has been implicated in obesity-related breast cancers (54). While clinical data relating to BMI or obesity status was not available for these samples, it is possible that this clustering could indicate obesity-related cases. Notably, the smaller subset of normal tissues corresponds to a significantly greater proportion of patients having received hormone therapy for breast cancer treatment (P < 0.05), indicating a potential association between such therapies and HERV activation. This is consistent with previous work in which HERV-K has been shown to be hormonally activated in breast cancers. There was no significant difference in other clinical variables between the two clusters, nor was there a difference in survival probabilities (P = 0.1), although this analysis may have been limited by the small sample size. As in prostate cancer, many of the differentially expressed genes intersect or are located in close proximity to differentially expressed genes (387 and 157, respectively), and several breast cancer-associated genes were found to be differentially expressed in these data, including BRCA1 and BRCA2. As in prostate cancer, several TET and APOBEC genes were also differentially expressed. One upregulated HERV (HARLEQUIN 17q21.31) intersects BRCA1, which was also upregulated in tumor samples (P < 0.05). In addition, NRN1—which is typically hypermethylated in breast tumors—is 3′ of downregulated element HERVH 6p25.1 (P < 0.01). In total, 1,094 of 1,547 differentially expressed HERVs intersected either a protein-coding gene or lncRNA.

In the TCGA-COAD dataset, 575 differentially expressed HERVs and 12,701 differentially expressed genes were identified. PAB and ERV1 groups were enriched among upregulated HERVs (q < 0.25). The separation of normal samples into two groups that was found in prostate and breast cancers was not seen in colon cancer, although this is possibly a result of the much lower sample size. Two tumor samples do separate along PC2 in HERV expression data, although this is not matched in hierarchical clustering or in PCA of gene expression data. Overall, differences in HERV expression between tumor and normal samples are less pronounced than in the other two datasets. Still, many differentially expressed HERVs intersect differentially expressed genes (119 HERVs) or are immediately 5′ or 3′ of a differentially expressed gene (128 HERVs). Many genes previously associated with colon cancer were identified as differentially expressed in these samples. As in the other cancer types, several differentially expressed HERVs intersect clinically relevant genes, and a majority (409/575) of HERVs intersect a protein-coding gene or lncRNA.

In total, 675 unique HERV loci were differentially expressed in more than one of the cancer types studied, suggesting that a significant amount of HERV expression in cancer is not disease or tissue specific. Moreover, the same HERV groups—MER4, HERVH, ERVL/HERVL, and HERVK—were most prevalent among dysregulated HERVs in all three cancer types, which is consistent with other studies identifying overlap in expressed HERV groups across TCGA cancer types (21). As the same groups were predominant among both up- and downregulated HERVs, collapsing HERV transcripts at the gene level may mask differential expression patterns. HERVK and HERVH were found to be two of the most frequently overexpressed HERVs in a previous study as well, where it was postulated that their overexpression was a result of being more recently integrated into the human genome and thus more functional (19). However, MER4 elements, which are older and thus no longer functional, were most prevalent. These results suggest that other HERV groups in addition to HERV-K should also be a focus of disease studies. Interestingly, the highest overlap of differentially expressed HERV loci is seen between prostate and breast cancers. These cancers commonly are associated with sex hormone dependence and mutation of BRCA1, suggesting that these two factors may also be linked to HERV expression. Furthermore, the subgroup ERV316A3 was repeatedly implicated in top differentially expressed HERVs in all three cancer types, motivating future study of these elements. Across all three datasets, many differentially expressed HERVs intersect genes relevant to disease phenotypes, and many are located in enhancer regions, motivating further studies into the functional effects of HERV expression in cancer tissues.

Demethylation is a key mechanism of HERV expression. Here, TET and APOBEC genes—which are important to DNA demethylation pathways—were differentially expressed in all three cancer types. Previous work has associated TET2 overexpression with increased HERV expression and antiviral immune response signatures (27), suggesting a link between demethylation and HERV activation in cancer. TET2 was downregulated in COAD, while TET3 was upregulated in BRCA and PRAD, with TET1 also being upregulated in PRAD. Thus, these results add to existing evidence suggesting changes to DNA demethylation activity contribute to HERV activation in tumor cells. Interestingly, both TET and APOBEC genes were both up- and downregulated in tumor samples, suggesting that epigenetic alterations in normal tissues are of importance, as well as tumor samples. Moreover, several genes associated with viral mimicry and host antiviral response – including MAVS, IRF7, and several interferon genes – were significantly overexpressed in tumor tissues (P < 0.05), which is consistent with previous work establishing demethylation of HERVs as a precursor to an immunogenic response (55). In addition, gene ontology enrichment analysis identified enrichment of several relevant biological pathways in dysregulated genes across cancer types, including epigenetic gene regulation and immune response. Therefore, gene expression data support the roles of both demethylation and viral mimicry pathways in contributing to HERV pathogenesis in cancer.

While further empirical studies will be required to validate the functional effects of HERV expression patterns identified here, significant empirical evidence exists for the immunogenic effects of HERV expression. In fact, HERVL/ERVL, HERVK, and ERV1—three HERV classes that were repeatedly identified in our analyses—were recently implicated in biochemical experiments linking HERV expression to interferon response in ovarian cancer cell lines (56, 57). Furthermore, HERV elements at several loci identified in the top significantly differentially expressed HERVs in our analyses have been identified in previous empirical studies, including the 21q22.3 locus (55). Finally, a recent study found evidence of HERV-induced transcriptional modulation of zinc-finger proteins in lung cancer cell lines, many of which were differentially regulated across all three cancer types studied here (21). Thus, previous empirical work supports the utility of our bioinformatics analyses in identifying targets of future validation studies.

The primary limitations of this study stem from low sample size. While paired tumor-normal tissue samples are necessary to identify dysregulated elements while controlling for heterogeneity among patients, normal solid tissue samples are difficult to obtain and thus are not commonly included in sequencing databases. In the TCGA-PRAD data in particular, survival data were not available, and thus the impact of HERV expression clustering on survival probability was not able to be determined. Also, only one patient in the TCGA-PRAD dataset had received hormone therapy, and so the effects of such treatment on HERV activity could not be evaluated for this dataset. In the TCGA-BRCA data, while survival data were available, the cluster of interest was small (n = 14), and so re-evaluation with a larger sample set would be beneficial. Moreover, all tumor samples included in this study were primary tumors. Additional samples from metastatic tumors would be useful to determine differential expression of HERVs throughout cancer progression. Finally, methylation data were not incorporated into this analysis as the methylation arrays used in the TCGA project primarily cover genic regions; however, further studies incorporating appropriate methylation data would be of great benefit toward uncovering causal relationships in HERV transcription pathways. While other studies have addressed HERV expression in cancer RNA-seq data, our study is unique in identifying locus-specific expression, as opposed to more general HERV groups or subgroups. In addition, the Telescope algorithm utilizes a greater proportion of reads by not requiring reads to be uniquely mapped via its Bayesian reassignment approach (28).

In conclusion, by utilizing RNA-seq data from matched samples, we identified differentially expressed HERVs in both tumor and normal tissue in prostate, breast, and colon cancers. Our work is novel in providing locus-specific characterization of these elements, which has previously been unattainable due to bioinformatics challenges. We found that HERV expression in normal tissue samples from prostate and breast cancers differentiates samples into two distinct groups, and that in breast cancer, this clustering was associated with hormone therapy. Furthermore, we found that a sizable number of HERV loci are differentially expressed across multiple cancer types, particularly between prostate and breast cancers, suggesting that the pathogenic role of HERVs in cancer may not be entirely disease-specific. By pairing HERV expression results with an analysis of host genes, we found that many DNA demethylation and antiviral response related genes were also differentially regulated in tumor tissues, providing evidence toward demethylation and viral mimicry as driving factors for HERV expression in cancers. Future work validating these findings using empirical data and determining causal relationships would be of great benefit, as such knowledge would elucidate regulatory mechanisms of disease-associated HERVs toward the development of novel cancer therapies.

M.C. Steiner reports grants from NIH during the conduct of the study. M.L. Bendall reports grants from NIH/NCI during the conduct of the study. K.B. Chiappinelli reports personal fees from Rome Therapeutics outside the submitted work. D.F. Nixon reports grants from NIH during the conduct of the study. K.A. Crandall reports grants from NIH during the conduct of the study. No disclosures were reported by the other authors.

M.C. Steiner: Conceptualization, resources, formal analysis, supervision, funding acquisition, investigation, visualization, writing–original draft, project administration, writing–review and editing. J.L. Marston: Funding acquisition, methodology, writing–review and editing. L.P. Iñiguez: Methodology, writing–review and editing. M.L. Bendall: Software, formal analysis, supervision, investigation, visualization, methodology, writing–original draft, writing–review and editing. K.B. Chiappinelli: Conceptualization, supervision, funding acquisition, methodology, project administration, writing–review and editing. D.F. Nixon: Conceptualization, supervision, funding acquisition, methodology, project administration, writing–review and editing. K.A. Crandall: Conceptualization, resources, software, supervision, funding acquisition, methodology, project administration, writing–review and editing.

This work was supported by the NIH, grant numbers CA206488 (to D.F. Nixon and K.A. Crandall), UL1TR000075 (to K.A. Crandall), and a Medical Scientist Training Program grant to the Weill Cornell-Rockefeller-Sloan Kettering Tri-Institutional MD-PhD Program T32GM007739 (to J.L. Marston). The work also benefited from support from a Research Innovation Award (02RIA032020) from the Milken Institute School of Public Health (to K.A. Crandall, K.B. Chiappinelli). The results shown here are based upon data generated by the TCGA Research Network: https://www.cancer.gov/tcga. We thank Miguel Branco and an anonymous reviewer for helpful comments to improve our manuscript.