Between May 2020 and June 2023, 20 patients diagnosed with malignant breast cancer who had sufficient tumor samples for whole-genome methylation sequence (GM-seq) were retrospectively enrolled in this study. The protocol was approved by the Ethics Committee of Zhongshan City People’s Hospital (K2023-198). Informed consent was obtained from all patients prior to sample collection and data use. ER/PR positivity in the immunohistochemistry testing were indicated with a cut-off of equal or higher than 1% positively stained cells. HER2 status was defined by an immunohistochemistry score of 3+ or presence of HER2 amplification (HER2/CEP17 ratio ≥ 2.0) demonstrated by fluorescence in situ hybridization (FISH) analysis. All patients received trastuzumab plus pertuzumab-based neoadjuvant therapy and subsequent radical surgery. pCR was defined as no evidence of the invasive cancer with or without intraductal components, including the lymph nodes, based on the pathological evaluation of the surgical breast specimen. A total of 28 perioperative tumor samples were obtained from 20 HER2-positive breast cancer patients treated with neoadjuvant therapy, including 12 pre-treatment samples from patients who achieved pCR, 8 pre-treatment samples and 8 matched post-treatment samples from 8 patients who did not (non-pCR).

Tissue DNA was extracted from FFPE by using QIAamp DNA MiniKit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. A total of 1 µg DNA was fragmented into 200–250 bp segments, using a Covaris S2 instrument (Woburn, MA, USA). DNA concentration was measured by Qubit™ dsDNA HS Assay Kit (ThermoFisher, Waltham, MA, USA).

We utilized the Hieff NGS Ultima Pro DNA Library Prep Kit for Illumina (Yeason, catalog number 12201ES96) for end repair, A-tailing, and adapter ligation. After magnetic bead purification, 5-methylcytosine (5mC) and 5-hydroxymethylcytosine (5hmC) are oxidized to 5-acylcytosine (5fC) or 5-carboxycytosine (5caC) by TET enzyme, and then used Pyridine borane treats 5fC or 5caC, reducing it to dihydrouracil (DHU) (16). Finally, PCR amplification was performed and barcodes were introduced to obtain a sequencing library. DNA libraries were assessed using Agilent Bioanalyzer. After this series of processes, the methylated cytosine was identified as thymine (T) by DNBSEQ-T7 (GenePlus-SuZhou, SuZhou, China) sequencing platform. Sequencing was performed with 150 bp paired-end reads, targeting >30 × coverage per sample. All libraries passed QC with >90% base quality score above Q30.

The analysis of DMRs between groups was conducted using Metilene (v0.2–8) (17), and only data with a CpG coverage depth greater than 5 were included. DMRs were required to contain a minimum of 10 CpG sites with ≤300 bp spacing between adjacent sites. Significant DMRs were identified based on two criteria: (a) absolute methylation difference > 0.1 between groups, and (b) q-value < 0.05. Functional annotation of DMR-associated genes was performed through GO enrichment analysis using the online website Metascape (https://metascape.org/gp/index.html). Methylation patterns were visualized via hierarchical clustering with pheatmap (v1.0.12), employing Euclidean distance metrics.

To ensure data reliability, only CpG sites with sequencing coverage ≥5× were retained for downstream genomic annotation. Gene promoter regions were defined as 1,000 bp upstream of transcription start sites (TSSs). Using BEDTools (v2.30.0), we extracted all qualifying CpG sites within these promoter regions. The promoter methylation level for each gene was then calculated as the mean β-value (methylation ratio = methylated reads/total reads) of all CpG sites within its promoter region. In the study, CpG sites were divided into three categories based on methylation levels: hyper (high methylation, beta value greater than 0.7), hypo (low methylation, beta value less than 0.3), and inter (moderate methylation, the rest of the range) (18).

All analyses were performed using R v4.2.0. Statistical comparisons were two-sided and adjusted using false discovery rate (FDR) where appropriate. Data visualization was conducted using ggplot2 and Complex Heatmap packages. Comparisons between two groups were based on two-sided Wilcoxon rank-sum (Mann–Whitney U) test, and comparisons among multiple groups were based on Kruskal-Wallis test.

20 female patients with stage II–III HER2-positive breast cancer were retrospectively enrolled in this study ( Table 1 ). The median age of the study cohort was 48 years (range 32–68 years). Most patients were HR-positive and HER2-positive (65%) and all patients were clinically lymph node positive. All of the involved patients received trastuzumab plus pertuzumab-based neoadjuvant therapy and radical surgery. 12 (12/20, 60%)patients achieved pCR, while 8 (8/20, 40%)patients were non-pCR after neoadjuvant therapy.

We collected 28 HER2-positive breast cancer samples from 20 patients who underwent neoadjuvant therapy, including 12 pre-treatment samples from patients who achieved pCR, 8 pre-treatment samples and 8 matched post-treatment samples from 8 patients who did not (non-pCR). In order to further characterize the molecular characteristics of HER2-positive breast cancer treatment response, whole-genome methylation sequencing was performed on all the above samples using sufficient tissues. For the CpG sites in each sample, they were classified based on their methylation levels, with sites with beta values greater than 0.7 defined as hyper, sites with beta values less than 0.3 defined as hypo, and the rest defined as inter. The average proportion of hyper sites in all samples was 0.39, the average proportion of hypo sites was 0.35, and the average proportion of inter sites was 0.26 ( Figure 1A ). In the pre-treatment samples, the non-pCR group had slightly more hypomethylation events than the pCR group. In the paired non-pCR samples, the post-treatment group had more hypermethylation events than the pre-treatment group ( Figure 1B ). At the same time, the differentially methylated sites could not only clearly distinguish pCR and non-pCR samples, but also distinguish non-pCR patients before and after treatment ( Figures 1C, D ), emphasizing the potential impact of methylation changes on the treatment response of HER2-positive breast cancer.

To further explore the methylation changes in the neoadjuvant therapeutic effect of HER2-positive breast cancer, we first performed differential methylation analysis in tumors before treatment. By comparing the tumors in the pCR and non-pCR groups, a total of 192 differentially methylated genes were identified. In non- pCR tumors, 61 hypomethylated and 131 hypermethylated gene promoters were identified separately. Hierarchical clustering analysis reveals distinct methylation profiles between the two groups, with non-pCR tumors exhibiting widespread hypermethylation compared to pCR tumors ( Figure 2A ). Correlation of methylation in differentially methylated regions showed that methylated genes in the same group tend to have co-hypermethylation or demethylation events ( Figure 2B ). The highly methylated genes in non-pCR tumors were enriched in multiple immune-related pathways, such as positive regulation of leukocyte activation, positive regulation of T cell activation, etc. ( Figure 2C ), which suggest that non-pCR tumors have immune disorders. Studies have shown that the HER2-L755S mutation can enhance HER2 autophosphorylation and affect downstream MAPK and PI3K/AKT/mTOR signaling pathways (19). ANKRD44 gene silencing activates NF-κB through the TAK1/AKT signaling pathway, and glycolysis levels increase in resistant cells (20). Interestingly, hypomethylated genes in non-pCR tumors were enriched in pathways such as phosphorylation (such as FGFR4, KIT, MOS, RET and ULK1) and canonical glycolysis ( Figure 2D ). This suggests that changes in gene promoter methylation may play a key role in the drug resistance mechanism that occurs with changes in different signaling pathways. Next, we focused our attention on oncogenes, which are genes that may be associated with the drug resistance mechanism. MOS, as a proto-oncogene, is also a serine/threonine kinase that activates the MAPK cascade by directly phosphorylating the MAP kinase activator MEK. Compared with pCR tumors, the level of methylation of the MOS gene in non-pCR tumors had a non-significant increasing trend ( Figure 2E ), suggesting that methylation variations in oncogenes may associated with the effect of treatment response. Another proto-oncogene that has attracted much attention is RET. Changes in the RET gene mainly include gene fusions and point mutations. Targeted therapies for RET gene changes mainly focus on RET kinase inhibitors, which mainly target RET gene fusions and certain point mutations, rather than methylation. However, changes in methylation cannot be ignored. In non-pCR tumors, there was a significantly lower methylation level, with a concentrated distribution and an average methylation level of only 0.17. In contrast, in pCR tumors, RET methylation levels were higher and relatively dispersed ( Figure 2F ). These phenomena may suggest that the correlation between treatment response and tumor methylation level.

Methylation changes within tumors after treatment still require further study. Therefore, we performed differential methylation analysis in paired samples before and after treatment in non-pCR patients. Compared with pre-treatment tumors, 284 significantly hypermethylated genes and 228 significantly hypomethylated genes were identified in tumors after treatment ( Figure 3A ). Similarly, methylated genes in the same group tend to undergo co-hypermethylation or demethylation events ( Figure 3B ). Among them, hypermethylated genes were involved in T cell activation, leukocyte activation, and other pathways involved in immune response ( Figure 3C ), while hypomethylated genes were involved in mesenchymal cell proliferation, positive regulation of MAPK cascade, and epithelial cell proliferation. ( Figure 3D ). These results demonstrated that treatment-resistant tumors might employ a biphasic epigenetic strategy: establishing immune privilege through hypermethylation of immunoregulatory genes while acquiring malignant potential via hypomethylation-activated mesenchymal and proliferative programs. To further understand the interactions between differentially methylated genes, we drew undirected interaction networks for hypermethylated and hypomethylated genes separately. It was worth noting that in the hypermethylated gene interaction network, we identified two key hub nodes, namely ERBB2 and PXDN ( Figure 3E ). The ERBB2 gene is a member of the epidermal growth factor receptor (HER) family and is related to a variety of cell signaling pathways that play important roles in cell proliferation, tumor formation, and apoptosis. In our study, no difference in ERBB2 methylation was found between pCR and non-pCR pre-treatment tumors, but increased promoter methylation of ERBB2 was found in non-pCR patients after treatment ( Figure 3F ). While increased promoter methylation may indicate gene silencing, this study does not directly measure ERBB2 expression, and further investigation is required. PXDN encodes a heme-containing peroxidase secreted into the extracellular matrix. It is involved in the formation of the extracellular matrix and may play a role in physiological and pathological fibrotic responses in the fibrotic kidney and can promote laminin assembly. Hypermethylation of Hub node genes may lead to decreased expression, thereby affecting downstream pathways. Hypermethylation of genes related to downstream pathways occurs at the same time, which further aggravates the drug resistance of tumors.

We hypothesized that genes exhibiting consistently increasing or decreasing promoter methylation levels across the three comparison groups (pCR, non-pCR before, and non-pCR after) would facilitate the identification of reliable resistance biomarkers. Therefore, we integrated the differentially methylated genes between pCR before treatment and non-pCR after treatment, and collected the RNA-sequencing profile of 285 HER2-positive breast cancer tumors receiving neoadjuvant therapy in GSE243375 dataset for verification. Finally, four genes with the same methylation difference trend were identified, namely KIT, LAD1, FAM110C, and DAPP1. As shown in Figure 4A , KIT exhibited progressively decreasing methylation levels in treatment-resistant cases, representing the only gene with this declining trend. In contrast, LAD1 (p=0.0013), FAM110C (p=0.0019) and DAPP1 (p=0.00036) all showed significant gradual methylation increases among the three groups. Notably, LAD1-encoded basement membrane anchoring protein showed reduced expression in poor responders ( Figure 4B , p=0.036). Similarly, DAPP1 - related to innate immunity - exhibited both progressively increased methylation and downregulated expression ( Figure 4B , p=0.03). These findings suggested that the methylation levels of genes related to epithelial-mesenchymal homeostasis and immune environment homeostasis may affect the treatment response of HER2-positive breast cancer.