Expression of Human Endogenous Retroviruses (HERVs) in Breast Cancer: A Transcriptomic Study

1. Introduction

Breast cancer remains a leading cause of cancer-related mortality among women globally [1]. In 2024, it is estimated that approximately 310.720 new cases of invasive breast cancer were diagnosed in women in the United States, along with about 56.500 cases of ductal carcinoma in situ (DCIS). Breast cancer continues to account for around 30% of all new cancer diagnoses in U.S. women. Mortality remains significant, with an expected 42.250 female deaths and 530 male deaths out of an estimated 2.790 new male cases. Despite these numbers, the five-year relative survival rate for breast cancer diagnosed at a localized stage exceeds 99%. This rate declines to 87% for regional-stage disease and drops further to 32% when diagnosed at a distant, metastatic stage. Encouragingly, breast cancer mortality has declined by 44% since 1989, resulting in over 517.900 averted deaths. However, concerning trends include a rising incidence among younger women under 50 years, increasing at an annual rate of 1.4% between 2012 and 2021 [2]. Significant racial disparities also persist, with Black women experiencing a 40% higher mortality rate compared to White women, largely due to later-stage diagnoses and unequal access to quality care [3].

Tumorgenesis is a multistep process, with oncogenic mutations in a normal cell conferring clonal advantage as the initial event. However, despite pervasive somatic mutations and clonal expansion in normal tissues, their transformation into cancer remains a rare event, indicating the presence of additional driver events for progression to an irreversible, highly heterogeneous, and invasive lesion [4]. Advances in molecular biology and genomics have revolutionized our understanding of tumor biology, yet many mechanisms contributing to breast cancer onset and progression remain poorly understood. Researchers are focusing on encompassing diverse aspects such as tumor stemness, intra-tumoral microbiota, and circadian rhythms [5]. Among the genomic elements receiving increasing attention are endogenous retroelements—particularly human endogenous retroviruses (HERVs) [6].

HERVs are sequences derived from ancient retroviral infections that have been integrated into the human genome over millions of years. Although most are defective and transcriptionally silent due to mutations and epigenetic repression, certain HERV families, especially human endogenous retrovirus type K [HERV-K (HML-2)], retain the ability to be transcribed under specific physiological and pathological conditions. In cancer biology, aberrant reactivation of HERVs has been reported in several malignancies, including melanoma, testicular cancer, and breast cancer [6,7,8].

The HERV-K (HML-2) subgroup is the most biologically active and has been implicated in oncogenesis through mechanisms such as insertional mutagenesis, induction of inflammation, and modulation of the immune microenvironment. The transcriptional activity of HERVs in tumor tissues may provide insights into novel oncogenic pathways and serve as a source of tumor-specific antigens [10,11]. The HERV-K (HML-2) is found to be expressed in various subtypes of breast cancer. This retrovirus was integrated into the primate genome approximately 55 million years ago, presenting a challenge for potential immunotherapy approaches for breast cancer. Recent studies have indicated that several endogenous retroviruses (ERVs) have been re-activated in tumors, with many exhibiting overexpression in cancerous tissues while showing low or undetectable levels in normal tissues [11,12].

This study focuses on evaluating the differential expression of HERVs in breast cancer versus normal breast tissue using publicly available RNA-Seq datasets. We also assessed the effect of treatment on HERV expression levels, with the ultimate goal of exploring HERVs as potential biomarkers and therapeutic targets in breast cancer.

2. Materials and Methods

2.1. Dataset Acquisition

RNA-Seq datasets were obtained from the NCBI Sequence Read Archive (SRA) [13]. All data has been securely downloaded into password-protected directories of the SRA. The SRA Toolkit was obtained through Linux to facilitate access to SRA data. This toolkit enables the creation of next-generation sequencing files in the specified format and cloud storage.

RNA-Seq data for breast ductal carcinoma (BAM files) along with their corresponding clinical information were sourced from the NCBI SRA. The paired-end FASTQ files for each sample were downloaded in SRA format and subsequently converted to FASTQ format using the fastq-dump command from the SRA Toolkit. Samples included 11 breast cancer tissues (various subtypes including ductal carcinoma and triple-negative) and 6 healthy breast tissues.

2.2. Bioinformatics Analysis

The paired-end reads in FASTQ format were aligned to the human reference genome, hg19 using bowtie2 with default settings The paired-end reads in FASTQ format were also aligned to the pseudogenome which was made using the 5’ and 3’ LTRs (long terminal repeat) coordinated of HERV-K(HML-2) as described by Subramanian et al. [14]. Mapping to this pseudogenome allowed for quantification of reads aligning to HERV-host junctions.

The resulting alignments were subsequently converted to the standard BAM file format utilizing samtools with its default configurations. The expression of HERV-K-human junctions was obtained from raw reads through bedtools multicov, also employing default settings. Following this, the mapped reads were quantified and normalized using the specified formula:

$${\displaystyle \frac{\mathit{m}\mathit{u}\mathit{l}\mathit{t}\mathit{i}\mathit{c}\mathit{o}\mathit{v}\mathit{p}\mathit{s}\mathit{e}\mathit{u}\mathit{d}\mathit{o}\mathit{g}\mathit{e}\mathit{n}\mathit{o}\mathit{m}\mathit{e}\mathit{r}\mathit{e}\mathit{a}\mathit{d}\mathit{s}}{\mathit{a}\mathit{l}\mathit{i}\mathit{g}\mathit{n}\mathit{m}\mathit{e}\mathit{n}\mathit{t}\mathit{r}\mathit{a}\mathit{t}\mathit{e}\mathrm{\%}\mathrm{*}\mathit{r}\mathit{e}\mathit{a}\mathit{d}\mathit{s}}}$$

A pseudogenome was constructed using the coordinates of LTRs from HERV-K(HML-2) to identify the expression of LTR-host junctions [15]. The LTR coordinates utilized for the development of the pseudogenome were detailed by Subramanian et al. The incorporation of HERVs into the human genome leads to areas that consist of both partial LTR sequences and partial human sequences, known as LTR-host junctions. Various tools were employed to analyze these datasets, including SRA-toolkit, Samtools, Bowtie, IGV, and BedTools, which is a software specifically designed to quantify the number of reads associated with each gene.

To incorporate the integration sites for the pseudogenome, partial human sequences were utilized, commencing and concluding at a defined distance of bases from the beginning and end of LTRs (the coordinates for the start and end of LTRs are provided in the attached Excel file). Consequently, a Bed file was generated by adjusting the length of reads by 66% to encompass junctions. The formula is outlined below:

chr(1,2,3, etc) LTR start-66% length LTR start+66% length

chr(1,2,3, etc) LTR end-66% length LTR end+66% length

for any LTR of HERV.

Command “bedtools getfasta” [16] was used with default settings, in order to create a fasta file including the LTR-host junction sequences that was used as the pseudogenome. As far as it concerns the alignment of pseudogenome, use the command:

“bowtie2 -t -x bowtie2/hg19 ncbi/public/sra/....fastq -S bowtie2/....sam”

Furthermore, downloading SAMtools was used for mapping the pseudogenome, via the following commands:

• ”wget https://github.com/samtools/samtools/releases/download/1.3.1/samtools-1.3.1.tar.bz2”
• “tar -jvxf samtools-1.3.1.tar.bz2”
• ”cd samtools-1.3.1”
• “make”

“.sam file” was converted to “.bam file” by means of the commands:

“samtools view -S -b bowtie2/....sam > bowtie2/....bam”

After that, files were sorted with samtools by means of the commands:

“samtools sort bowtie2/....bam -o bowtie2/..._sorted.bam”

“samtools index bowtie2/..._sorted.bam” [17].

2.3. Statistical Analysis

Data analysis was performed in IBM SPSS 24. Chromosomes were represented using frequencies. Normalized coverage of integration sites was calculated using BEDTools and custom scripts. Statistical significance between groups was assessed using independent samples t-tests and Mann-Whitney U tests, with significance set at p < 0.05.

3. Results

3.1. RNASeq Datasets

Transcriptome data were examined from 11 patients with invasive ductal carcinoma breast cancer, as detailed in Table 1, alongside 6 healthy individuals, presented in Table 2.

In Table 3 and Table 4, the raw reads and overall alignment rate from the alignment of any RNA-seq to Human genome (hg19) are presented for breast cancer patients and healty tissue respectively. Pseudogenome reads came from alignment of any RNA sequences to Pseudogenome. Normalized coverage of integration sites:

$${\displaystyle \frac{multicovpseudogenomereads}{alignmentrate\mathrm{\%}\mathrm{*}reads}}$$

3.2. Elevated HERV Expression in Breast Cancer

Mean normalized coverage of HERV integration sites was significantly higher in breast cancer tissues (M = 0.000013) than in healthy tissues (M = 0.000001, p = 0.009) (Table 5).

3.3. Effect of Treatment

Patients who had not received treatment exhibited significantly higher HERV expression (mean rank = 9.00) compared to treated patients (mean rank = 4.29, p = 0.023) (Table 6).

3.4. Visualization

The visualization of files was conducted using IGV, a software application that enables the observation of reads on chromosomes [17,18]. The associated IGV web application enables users to upload FASTQ files and receive visualized outcomes. IGV genome browser confirmed increased read accumulation at specific HERV loci in cancer samples versus healthy controls (Figure 1, Figure 2, Figure 3 and Figure 4).

4. Discussion

The present analysis investigated the expression of HERVs in breast cancer tissue compared to normal breast tissue. Furthermore, it explored the variations in normalized coverage of integration sites between breast cancer patients and healthy controls, as well as the differences between treated and untreated breast cancer cases.

The results provide compelling evidence that HERVs, specifically HERV-K (HML-2), are significantly overexpressed in breast cancer tissues compared to normal breast tissues. This differential expression pattern supports a growing body of literature suggesting that HERVs play a role in the etiology and progression of cancer [19,20,21].

Several mechanisms may explain the upregulation of HERVs in malignancies. Hypomethylation of HERV sequences and surrounding genomic regions, a common epigenetic alteration in cancer, may reactivate normally silenced proviral elements. Additionally, the tumor microenvironment, characterized by oxidative stress, cytokine production, and immune dysregulation, can further promote HERV expression [22,23]. Specifically in breast cancer, the envelope protein of HERV-K(HML-2) serves as the most effective diagnostic tool for breast cancer, as it is found in breast cancer tissues while being absent or present at significantly lower levels in normal breast tissue [24]. Also, the presence of HERV-K (HML-2) antibodies and mRNA serves as an indicator of early-stage breast cancer, with a subsequent rise suggesting the occurrence of metastasis [20,25]

The observed reduction in HERV expression among patients who received treatment is particularly noteworthy. It suggests that chemotherapy or other therapeutic regimens may suppress HERV activity, potentially by restoring epigenetic silencing mechanisms or by targeting proliferative cancer cells with elevated HERV transcription [19,26,27]. This raises important questions regarding the therapeutic targeting of HERVs, both directly and indirectly [28].

While this study offers valuable insights, limitations include the modest sample size and the reliance on public datasets without experimental validation. Further research is necessary to elucidate the functional consequences of HERV expression and its interaction with oncogenic signaling pathways.

5. Implications for Practice and Future Directions

The findings of this study have several clinical and research implications.

In clinical practice, HERV-K expression could serve as a non-invasive biomarker for early breast cancer detection [24]. The observed decrease in HERV expression following treatment suggests potential utility in tracking therapeutic efficacy. Furthermore, stratifying patients by HERV expression profiles may aid in developing tailored treatment strategies and enhancing personalized medicine approaches. The field of personalized medicine is experiencing a significant transformation due to the incorporation of multi-omics data, which primarily includes genomics, transcriptomics, proteomics, and metabolomics [29].

From a research perspective, laboratory-based validation is essential to confirm the functional role of HERVs in tumor growth and immune modulation. The immunogenic nature of HERV-derived proteins also presents opportunities for the development of novel immunotherapies, including cancer vaccines and CAR-T cell therapy [30]. Lastly, longitudinal studies monitoring HERV expression over time could provide valuable insights into disease progression, recurrence, and metastasis.

6. Conclusions

This transcriptomic analysis demonstrates significant overexpression of HERVs in breast cancer tissues, particularly in untreated patients. HERVs expression profiling may aid in diagnosis, prognosis, and treatment response monitoring in breast cancer care. With further validation, HERVs could become valuable biomarkers and therapeutic targets in oncology.