This study integrated in silico transcriptomic discovery with clinical cohort validation. Public RNA-seq data from dataset GSE8427520 deposited in the Gene Expression Omnibus (GEO)21 were analyzed to prioritize dysregulated HERV-K (HML-2)-associated env signals in breast epithelial transformation models. Subsequently, family-level env expression was quantified via qRT-PCR in peripheral blood samples from BC patients and healthy controls. The clinical diagnostic validation component was conducted and reported in strict accordance with the Standards for Reporting Diagnostic Accuracy (STARD 2015) guidelines (Supplementary File S1).

Normalized RNA-seq data from GSE84275 were retrieved from GEO and evaluated across one non-transformed human mammary epithelial model (HME) and three transformed models (HMLE-Ras, HMLE-Her2, and HCC1954)20. HERV-K (HML-2)-associated env signals were identified based on repetitive element annotations and genomic coordinates overlapping known HERV-K loci. Relative expression differences between transformed and non-transformed models were evaluated descriptively using normalized expression values. Given the exploratory nature and limited sample size of this public dataset, formal hypothesis testing and multiple-testing corrections were not performed; the dataset was utilized primarily for candidate prioritization. The HERV-K (HML-2)-associated env signal overlapping chromosome 4p16.1 (chr4:9,123,514–9,133,075) demonstrated consistent upregulation across transformed models and was selected for clinical validation. Genomic coordinates were mapped to the GRCh38/hg38 reference assembly via the UCSC Genome Browser22,23.

The candidate HERV-K (HML-2) env-associated sequence was aligned against the human reference genome using BLASTn (NCBI) to verify sequence similarity and genomic overlap. Annotated genes situated within ±500 kb of the overlapping genomic region were extracted from the UCSC Table Browser using GENCODE v46 annotations24. Given the multicopy nature of HERV-K (HML-2) proviruses, this analysis was intended to confirm family-level sequence identity and genomic context rather than definitive locus-specific transcriptional assignment25.

Female patients with primary breast cancer were prospectively recruited from Hue University of Medicine and Pharmacy Hospital between January 2025 and January 2026. Diagnoses were histopathologically confirmed, and anatomical tumor staging was classified according to the American Joint Committee on Cancer (AJCC) TNM staging manual (8th edition)26. Early-stage disease was defined as stages I to IIA, and advanced-stage disease was defined as stages IIB to IV27.

Inclusion criteria were: (1) female sex aged ≥30 years; (2) histopathologically confirmed primary breast carcinoma; and (3) no prior history of systemic chemotherapy, radiotherapy, or endocrine therapy before blood collection. Exclusion criteria comprised: (1) previous history of other malignancies; (2) autoimmune or systemic inflammatory disorders; (3) acute infectious episodes within two weeks prior to sampling; and (4) pregnancy or lactation. Age-matched healthy controls were recruited from women undergoing routine annual health check-ups at the same institution during the corresponding timeframe. Control individuals exhibited no history of malignancy, autoimmune disease, or active infection. All participants provided written informed consent. The protocol was approved by the Institutional Ethics Committee of Hue University of Medicine and Pharmacy, Hue University (Approval No. H2025/426) and conducted in accordance with the Declaration of Helsinki (1975).

Baseline demographic characteristics and laboratory parameters, including complete blood counts and fasting blood glucose levels, were collected from electronic medical records. In BC patients, serum concentrations of cancer antigen 15-3 (CA15-3), carcinoembryonic antigen (CEA), and C-reactive protein (CRP) were determined through routine hospital laboratory assays. Fasting peripheral venous blood (4 mL) was collected in the morning: 2 mL in EDTA tubes for leukocyte isolation and RNA extraction, and 2 mL in non-anticoagulated tubes for biochemical testing. Samples were processed within 2 hours of venipuncture. Leukocytes were isolated using red blood cell (RBC) lysis buffer (Solarbio, Beijing, China) and stored at −80 °C until analysis.

Total RNA was extracted from isolated leukocytes using the AxyPrep Total RNA Miniprep Kit (Axygen, China). RNA concentration and optical purity (A260/A280 ratio) were assessed spectrophotometrically. Reverse transcription was performed using the HiScript® III RT SuperMix kit (Vazyme, China), which incorporates an initial genomic DNA elimination step. Quantitative real-time PCR (qRT-PCR) was executed using SYBR Green chemistry on a StepOne™ Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). Reactions were set up in a 20 μL volume containing 2× SYBR Green Master Mix, 0.4 μM forward and reverse primers, cDNA template, and nuclease-free water. All experimental procedures adhered strictly to the Minimum Information for Publication of Quantitative Real-Time PCR Experiments (MIQE) guidelines.

Target gene specificity was confirmed by post-amplification melting curve analysis, demonstrating single, discrete dissociation peaks (Supplementary Figure S1). Reactions were run in technical duplicates with no-template controls (NTC) and no-reverse-transcriptase (no-RT) controls in each batch. Relative HERV-K (HML-2)-associated env expression was calculated using the comparative 2^−ΔΔCt method, with glyceraldehyde 3-phosphate dehydrogenase (GAPDH) serving as the internal reference gene28. A single pooled cDNA calibrator derived from five BC patients was included across all runs to ensure batch-to-batch comparability: ΔCt = Ct(env) − Ct(GAPDH), and ΔΔCt = ΔCt(sample) − ΔCt(calibrator). Validated primers targeting conserved HERV-K (HML-2) env family sequences were utilized29 (amplicon size: 166 bp; Table 1).

Statistical analyses were carried out using SPSS 25.0 (IBM Corp., Armonk, NY, USA) and GraphPad Prism 9.0 (GraphPad Software, San Diego, CA, USA). Continuous variables were evaluated for distributional normality using the Shapiro–Wilk test and are presented as mean ± standard deviation (SD) for normally distributed data or median (interquartile range, IQR: 25th–75th percentiles) for non-normally distributed data. The sample size was exploratory and not predetermined.

Differences between two groups were assessed using the independent Student's t-test or the Mann–Whitney U test, as appropriate. Differences among three or more groups were evaluated using the Kruskal–Wallis test followed by Dunn's post hoc multiple comparisons test. Binary logistic regression analysis was conducted to identify factors independently associated with breast cancer. Multicollinearity among independent variables was assessed using variance inflation factor (VIF) and tolerance metrics. Receiver operating characteristic (ROC) curve analysis was performed to determine diagnostic accuracy, reporting the area under the curve (AUC) and 95% confidence intervals (CI). All statistical tests were two-tailed, and P < 0.05 was defined as statistically significant.

Exploratory screening of normalized RNA-seq profiles from GSE84275 demonstrated marked transcriptional elevation of an HERV-K (HML-2) env-associated signal overlapping the genomic interval chr4:9,123,514–9,133,075 in transformed breast epithelial cell lines (HMLE-Ras, HMLE-Her2, and HCC1954) relative to non-transformed mammary epithelial cells (HME). These transcriptomic findings served as a hypothesis-generating basis for clinical evaluation in patient-derived peripheral blood samples.

BLASTn alignment of the candidate sequence revealed 100% identity (7,836/7,836 bp, 0 gaps, E-value = 0.0) against the human reference provirus HERV-K (HML-2) sequence (JN675026.1; Supplementary Figure S2), corresponding to the HML-2_4p16.1a locus at chromosome 4p16.1 (chr4:9,123,514–9,133,075)20. Because HERV-K elements exhibit high sequence paralogy, this in silico confirmation established family-level identity and genomic overlap without inferring exclusive single-locus transcription.

Annotation via the UCSC Genome Browser (GRCh38/hg38) demonstrated that the candidate locus at 4p16.1 possesses a canonical full-length proviral architecture flanked by long terminal repeat (LTR5) sequences (Figure 1), known to harbor potent promoter and enhancer activities30. In addition, interrogation of the ±500 kb flanking region identified several pseudogenes and a microRNA gene (MIR548L2) (Supplementary Table S1).

A total of 101 female participants, comprising 56 BC patients and 45 age-matched healthy controls, were enrolled. Mean age did not differ significantly between BC patients and controls (51.50 vs. 47.00 years, P = 0.1632), confirming baseline comparability. BC patients exhibited significantly lower red blood cell counts (4.254 ± 0.4212 vs. 4.494 ± 0.2522 × 10¹²/L, P = 0.0011), lower hemoglobin levels (120.9 ± 12.06 vs. 133.6 ± 7.405 g/L, P < 0.0001), and a reduced lymphocyte percentage (28.23 ± 8.498% vs. 33.26 ± 6.312%, P = 0.0013). Conversely, platelet counts were moderately elevated in BC patients (265.5 vs. 235.0 × 10⁹/L, P = 0.0347). Total leukocyte count and fasting blood glucose did not differ between cohorts (all P > 0.05; Table 2, Supplementary Table S2).

qRT-PCR quantification revealed significantly higher peripheral blood leukocyte HERV-K (HML-2) env expression in BC patients than in healthy controls (median [IQR]: 0.9111 [0.5364–1.312] vs. 0.5284 [0.2848–0.7876], P = 0.0003; Figure 2). The Hodges–Lehmann median difference between cohorts was 0.3173.

Univariate logistic regression revealed that elevated HERV-K (HML-2) env expression (OR = 5.28, 95% CI: 1.90–14.67, P = 0.001), lower lymphocyte percentage (OR = 0.92, 95% CI: 0.86–0.97, P = 0.003), and lower hemoglobin level (OR = 0.87, 95% CI: 0.82–0.92, P < 0.001) were significantly associated with BC. In multivariable logistic regression, elevated env expression remained independently associated with BC (adjusted OR = 12.66, 95% CI: 3.14–51.02, P < 0.001). Patient age also demonstrated an independent association (adjusted OR = 1.05, 95% CI: 1.00–1.10, P = 0.044), while hemoglobin level maintained an inverse association (adjusted OR = 0.85, 95% CI: 0.79–0.91, P < 0.001). Lymphocyte percentage was not significant after adjustment (P = 0.179; Table 3). Multicollinearity diagnostics confirmed acceptable tolerance (0.832–0.996) and VIF values (1.004–1.203; Supplementary Table S3).

ROC curve analysis showed that peripheral HERV-K (HML-2) env expression effectively discriminated BC patients from healthy controls, achieving an AUC of 0.7075 (95% CI: 0.6051–0.8100, P = 0.0005; Figure 3), reflecting moderate diagnostic discriminatory accuracy.

Stratification of BC patients by disease stage (early-stage [I–IIA], n = 22; advanced-stage [IIB–IV], n = 34; healthy controls, n = 45) demonstrated significant stage-dependent heterogeneity (Kruskal–Wallis H = 21.54, P < 0.0001; Figure 4). Dunn's post hoc analysis indicated that env expression was significantly higher in early-stage BC patients than in healthy controls (adjusted P < 0.0001) and advanced-stage patients (adjusted P = 0.0166). In contrast, advanced-stage patients did not differ significantly from healthy controls (adjusted P = 0.1481). Median expression was highest in early-stage BC (1.087 [0.8586–1.448]), intermediate in advanced-stage BC (0.7196 [0.4563–1.104]), and lowest in controls (0.5284 [0.2848–0.7876]). Conventional markers (CA15-3, CEA, and CRP) did not differ significantly between stages (all P > 0.05; Supplementary Table S4). Clinicopathological characteristics and subtype distributions are summarized in Supplementary Table S5.