Classification of breast cancer cell lines into three distinct subgroups based on chromatin accessibility patterns

Using chromatin accessibility analysis, we investigated the epigenetic landscape in 23 breast cancer cell lines, including 2 ER+/HER2−, 3 ER+/HER2+, 3 HER2+, and 15 TNBC lines (Supplementary Table 1). All cell lines exceeded a predefined threshold for the transcription start site (TSS) enrichment score (≥5) and exhibited a unique fragmentation size distribution pattern with nucleosomal periodicity (Supplementary Fig. 1). Through peak calling, we identified 140,246 reproducible cis-regulatory elements (CREs) (Fig. 1a). A large proportion of these CREs were located in distal regions, including distal intergenic and intronic regions, whereas promoter elements constituted only 21.7% of the total CREs, consistent with previous findings [6, 7].

Fig. 1: Classification of breast cancer cell lines based on chromatin accessibility profiles.

a Genomic features of 140,246 merged reproducible peak sets. UTR: untranslated region. b Heatmap of Pearson’s correlations for ATAC-seq signals with all reproducible peaks. c Principal component analysis using the ATAC-seq profiles of cell lines. d Heatmap showing the chromatin accessibility of the top 5000 variable peaks. Annotations above the heatmap represent the receptor statuses of cell lines.

Based on CRE accessibility, correlation analysis among cell lines enabled the identification of three distinct groups: the receptor-positive group (Group-P), encompassing ER+ and/or HER2+ lines, including T47D, MCF7, MDA-MB-361, BT474, YMB1, MDA-MB-453, and SKBR3; the basal group (Group-B), primarily comprising basal-like TNBC lines, such as HCC1937 (BL1), MDA-MB-468 (BL1), HCC1806 (BL2), HCC38 (BL1), HCC1143 (BL1), and HCC70 (BL2) [8]; and the mesenchymal group (Group-M), including both mesenchymal and mesenchymal stem-like lines, such as BT549 (M), MDA-MB-436 (MSL), MDA-MB-231 (MSL), Hs578T (MSL), and MDA-MB-157 (MSL) [8] (Fig. 1b). Principal component analysis confirmed the segregation of these chromatin accessibility groupings (Fig. 1c). The cell line DU4475 displayed a unique chromatin accessibility pattern, likely due to distinct biological properties arising from its cutaneous metastatic nodule origin [15]. Notably, the accessibility pattern of the top 5000 most variable CREs indicated that Group-B exhibited high accessibility for CREs specific to both Group-P and Group-M (Fig. 1d). This suggests that Group-B has intermediate features of chromatin accessibility compared with the other groups. In total, 23 breast cancer cell lines were categorized into the three aforementioned groups based on chromatin accessibility.

Motif enrichment of FOXA1 and GRHL2 associated with chromatin accessibility–based groups of breast cancer cell lines

To gain insight into the three distinct subgroups, we analyzed enrichment of TF binding motifs within chromatin accessible regions across cell lines. Clear differences in motif enrichment patterns were observed between Group-P and Group-M (Fig. 2a). In particular, motifs belonging to the forkhead TF family, E2A, and GRHL2 were highly enriched in Group-P but less enriched in Group-M. In contrast, TEAD, RUNX, and AP-1 family motifs exhibited greater enrichment in Group-M, underscoring their crucial role in the mesenchymal state of breast cancer [16,17,18]. Notably, Group-B showed low enrichment of the forkhead motifs but displayed high enrichment of the GRHL2 motif (Fig. 2a).

Fig. 2: Enrichment of FOXA1 and GRHL2 binding motifs and chromatin accessibility in their coding regions.

a Heatmap representing ChromVAR motif scores for the top 30 most variable motifs. Annotations above the heatmap indicate chromatin accessibility groups and the receptor statuses of cell lines. Boxplots showing motif scores across chromatin accessibility groups for ER (b), FOXA1 (c), and GRHL2 (d). P-values, calculated via Student’s t-test, are shown. Genome track view of the ATAC-seq profiles of cell lines at loci around FOXA1 (e) and GRHL2 (f). MSL mesenchymal stem-like, M mesenchymal, Un unclassified, IM immunomodulatory, BL1 basal-like1, BL2 basal-like2. g Bar plots showing the relative expression of FOXA1 (blue) and GRHL2 (red) compared to MDA-MB-231. Error bars represent standard deviation.

We subsequently examined the luminal-lineage TFs ER and FOXA1. Enrichment of the ER-responsive element exhibited low variation among the three groups (Fig. 2b). In contrast, FOXA1 enrichment clearly distinguished Group-P from the other two groups (Fig. 2c), highlighting FOXA1’s importance as a TF associated with receptor-positive tumor lineages. Interestingly, GRHL2 motif enrichment remained consistently high in both Group-P and Group-B but was significantly lower in Group-M (Fig. 2d). These results suggest that the coexistence of FOXA1 and GRHL2 is unique to Group-P, whereas Group-B is characterized by the absence of FOXA1 and the presence of GRHL2. In contrast, Group-M is characterized by the absence of both FOXA1 and GRHL2.

To investigate the transcriptional activity of FOXA1 and GRHL2, we evaluated the accessibility of genomic loci near these two genes and analyzed their expression. Group-P cell lines showed substantial chromatin accessibility not only around the FOXA1 TSS but also in regions upstream and downstream of FOXA1, indicating the presence of enhancer elements (Fig. 2e). In contrast, Group-B cell lines exhibited moderate accessibility at the FOXA1 TSS but lacked accessibility in the enhancer regions observed in Group-P (Fig. 2e). The GRHL2 TSS exhibited high accessibility in Group-P and Group-B but lower accessibility in Group-M (Fig. 2f). Consistent with these observations, the mRNA expression levels of both FOXA1 and GRHL2 are high in Group-P cell lines (Fig. 2g). However, in Group-B cell lines, FOXA1 expression is lower than in Group-P, while GRHL2 expression remains at the same level as in Group-P lines. In Group-M, both genes show very low expression. These findings are consistent with the motif enrichment analysis and underscore the relationship between FOXA1 and GRHL2 activities and chromatin accessibility patterns in breast cancer cell lines.

Functional annotation of group-specific CREs reveals the distinct epigenetic landscape associated with each group’s unique properties

To explore epigenetic distinctions among groups, we conducted a differential accessibility analysis, resulting in the identification of six CRE sets: Group-P–specific (N = 8650), Group-B–specific (N = 1552), Group-M–specific (N = 7660), Group-P/B–shared (N = 8498), Group-B/M–shared (N = 12,505), and Group-M/P–shared (N = 490) (Fig. 3a; Supplementary Table 2). Consistent with prior findings indicating an intermediate epigenetic state in Group-B compared with the other groups (Fig. 1d), we observed a substantial number of accessible CREs shared between Group-P and Group-B, whereas the number of CREs shared between Group-P and Group-M was comparatively lower. Motif enrichment analysis revealed significant enrichment of FOXA1 and GRHL2 motifs in Group-P–specific and Group-P/B–shared CREs (Fig. 3b; Supplementary Table 3). Conversely, GRHL2 motifs were significantly enriched in Group-B–specific CREs, whereas FOXA1 motifs were absent.

Fig. 3: Group-specific or shared CREs.

a Heatmap showing the chromatin accessibility of group-specific or shared CREs. Annotations above the heatmap indicate chromatin accessibility groups and the receptor statuses of cell lines. b Annotation of each CRE set and significantly enriched motifs. c Bar plots of GO enrichment obtained using GREAT analysis of each CRE set.

We subsequently performed GREAT [19] Gene Ontology (GO) analyses for each set of CREs. We found that Group-P–specific CREs were associated with gland morphogenesis (Fig. 3c; Supplementary Table 4). Conversely, Group-M–specific CREs exhibited associations with processes related to extracellular matrix organization and epithelial cell migration, indicating their utility in characterizing the mesenchymal traits of Group-M cell lines. Notably, both Group-B–specific and Group-P/B–shared CREs showed associations with skin development, epidermal development, and keratinocyte differentiation. These findings suggest that group-specific CREs have a marked impact on establishing cellular identity. Collectively, these findings emphasize the critical role of FOXA1 in maintaining the epigenetic state of luminal or receptor-positive cells. Although both Group-B and Group-M are TNBC cell lines, they exhibit unique epigenetic profiles, with Group-B’s epigenetic state being notably influenced by GRHL2.

FOXA1 and GRHL2 play a pivotal role in regulating the accessibility of Group-P–specific, Group-B–specific, and Group-P/B–shared CREs

Although motif enrichment analysis offers valuable insight into the potential activity of TFs, it cannot definitively establish their roles. Therefore, we conducted knockdown experiments targeting FOXA1 (FOXA1-KD) and GRHL2 (GRHL2-KD) in T47D (Group-P) and HCC38 (Group-B) to determine their involvement in the chromatin accessibility of Group-P–specific, Group-B–specific, and Group-P/B–shared CREs. We examined the expression analysis after knockdown, which showed the successful knockdown of both genes in these cell lines (Supplementary Fig. 2a–d).

We first knocked down FOXA1 and GRHL2 in T47D followed by ATAC-seq analysis. FOXA1-KD led to reduced accessibility to Group-P–specific CREs (Fig. 4a, b) and, surprisingly, increased accessibility to Group-B–specific CREs (Fig. 4c). In contrast, the accessibility of Group-P/B–shared CREs remained unaffected (Fig. 4d). GRHL2-KD caused a slight decrease in the accessibility of Group-P/B–shared CREs, with no effects on either Group-P–specific and Group-B–specific CREs (Fig. 4a–d). These findings suggest that FOXA1 maintains the accessibility of Group-P–specific CREs while suppressing the accessibility of Group-B–specific CREs, with no significant effect on the accessibility of Group-P/B–shared CREs. Furthermore, GRHL2 is partially involved in regulating Group-P/B–shared CREs in the Group-P cell line.

Fig. 4: Effects of FOXA1 or GRHL2 knockdown on the chromatin accessibility of Group-P–specific, Group-B–specific, and Group-P/B–shared CREs.

a Heatmap showing normalized ATAC-seq signals in Group-P–specific, Group-B–specific, and Group-P/B–shared CREs of control, FOXA1 or GRHL2 knockdown in T47D (Group-P) cells. Normalized read count profiles in Group-P–specific (b), Group-B–specific (c), and Group-P/B–shared (d) CREs of control, FOXA1 or GRHL2 knockdown in T47D (Group-P) cells. e Heatmap showing normalized ATAC-seq signals in Group-P–specific, Group-B–specific, and Group-P/B–shared CREs of control, FOXA1 or GRHL2 knockdown in HCC38 (Group-B) cells. Normalized read count profiles in Group-P–specific (f), Group-B–specific (g), and Group-P/B–shared (h) CREs of control, FOXA1 or GRHL2 knockdown in HCC38 (Group-B) cells.

Subsequently, FOXA1- and GRHL2-KD was conducted in HCC38, followed by ATAC-seq analysis. For Group-P–specific CREs, we observed low accessibility of these regions in HCC38 under three conditions, specifically a low number of reads per million mapped reads. (Fig. 4e, f). For Group-B–specific CREs, a slight decrease in accessibility was observed with GRHL2-KD, while FOXA1-KD showed no change (Fig. 4g). Furthermore, GRHL2-KD resulted in a reduction in Group-P/B–shared CRE accessibility, whereas FOXA1-KD did not produce the same effect (Fig. 4h). These results suggest that GRHL2 plays a critical role in preserving the shared epigenetic signature between Group-P and Group-B cell lines.

FOXA1-KD reduces cell cycle in T47D cells, while GRHL2-KD activated a mesenchymal gene expression program in HCC38 cells

We investigated the effects of chromatin accessibility changes, following FOXA1 and GRHL2 knockdown, on gene expression regulation and phenotypes. RNA-seq analysis revealed 407 upregulated and 1029 downregulated genes in T47D following FOXA1-KD, and 183 upregulated and 146 downregulated genes following GRHL2-KD (Fig. 5a, b; Supplementary Table 5). In HCC38, GRHL2-KD resulted in 203 upregulated and 136 downregulated genes (Fig. 5c; Supplementary Table 5). Gene enrichment analysis revealed that genes linked to the cell cycle were downregulated by FOXA1-KD in T47D (Fig. 5d), highlighting the critical role of FOXA1 in cell proliferation in luminal breast cancer [20]. In T47D cells, GRHL2-KD transcriptionally upregulated genes associated with apoptosis. Furthermore, genes upregulated by GRHL2-KD in both T47D and HCC38 were associated with epithelial-mesenchymal transition (EMT) (Fig. 5d). Notably, GRHL2-KD in HCC38 upregulated well-known EMT markers such as VIM and VCAN, as well as a key EMT regulator, SNAI2 [21](Supplementary Table 5). These findings indicate that GRHL2 regulates gene expression, distinguishing between basal and mesenchymal lineages.

Fig. 5: Transcription and cell proliferation changes of FOXA1 or GRHL2 knockdown.

MA plots showing differentially expressed genes between the control and FOXA1 knockdown in T47D (a), GRHL2 knockdown in T47D (b), and GRHL2 knockdown in HCC38 (c). d Gene enrichment analysis for each set of differentially expressed genes. No enriched terms were found for downregulated genes by siGRHL2 in T47D and in HCC38. Bar plots showing cell viability of T47D (e) and HCC38 (f) after siRNA transfection for 48, 72, 96, 120 h. Error bars represent standard deviations. P-values, calculated via Student’s t test, are shown.

The cell proliferation assay showed that FOXA1-KD decreased cell proliferation in T47D cells, consistent with the decreased expression of cell cycle genes caused by FOXA1-KD (Fig. 5e). Similarly, GRHL2-KD resulted in decreased cell proliferation in T47D. In contrast, GRHL2-KD increased cell proliferation in HCC38 cells (Fig. 5f), suggesting a different role for GRHL2 in Group-P and Group-B cell lines.

Finally, we conducted survival analysis using The Cancer Genome Atlas invasive breast cancer (TCGA-BRCA) data. We found that elevated FOXA1 expression may be associated with poorer prognosis in luminal A tumors (P = 0.067, Log-rank test; Supplementary Fig. 3a, b), but no such correlation was observed in basal tumors. Furthermore, elevated levels of GRHL2 expression were significantly associated with a poorer prognosis in luminal A tumors (P = 0.038, Log-rank test; Supplementary Fig. 3c, d). However, this correlation was not statistically significant in basal tumors. Taken together, these findings suggest that FOXA1 plays a critical role in tumor proliferation in Group-P cells. Conversely, GRHL2 possesses a distinct function in maintaining cell viability in Group-P and sustaining the basal program in Group-B cells through the regulation of chromatin accessibility and gene expression.