Rank signaling drives basal cell-lineage infidelity leading to mammary tumorigenesis

Domain: chemistry
Published: 2026-02-27
Editorial outcome: major_revision

Abstract
Rank signaling regulates mammary gland development and epithelial differentiation. While Rank is expressed in both basal and luminal cells, its basal-specific role is unclear. Here, using inducible basal-specific Rank expression and lineage tracing, we show that Rank signaling regulates basal cell identity in postnatal mammary glands. Increased basal Rank activity disrupts basal and luminal identities, causing aberrant luminal-like differentiation, lactation defects, and premalignant lesions composed of hybrid basal-derived cells that progress to basal and luminal adenocarcinomas. Conversely, Rank loss reduces tumor formation and also impairs cell identity. Mechanistically, proteomic, transcriptomic, and chromatin analyses reveal that Rank activation drives epigenetic remodeling, leading to basal identity loss and tumor initiation. A basal Rank gene signature correlates with ductal carcinoma in situ recurrence, as well as poor outcomes in luminal breast cancers. Thus, basal Rank-driven lineage infidelity promotes pre-invasive lesions and transition to invasive breast cancer in females.

Introduction
The mammary gland (MG) is a specialized organ that develops postnatally and experiences dynamic morphological changes, from an initial expansion during puberty to subsequent cycles of growth/regression led by oestrous and pregnancy/involution cycles. The MG is organized as an epithelial bilayer comprising luminal cells (LCs) and basal cells (BCs) surrounding the LC layer. LCs are either hormone-sensing cells, which express hormone receptors (Hrs), or luminal progenitors (LPs), devoid of Hrs. Embryonic MG development relies on bipotent progenitors located in the basal compartment, able to give rise to LCs and BCs (reviewed in 1). Upon birth, each epithelial compartment is mostly maintained by lineage-restricted progenitors (reviewed in 1). However, postnatally BCs retain the competence for bipotency, triggered under specific conditions 2,3.

Rank (receptor activator of nuclear factor kappa-B) and its ligand, Rankl, members of the TNF (receptor) superfamily, play pivotal roles in MG development 4,5, being the role of Rank signaling in LCs extensively studied 3-8. Rankl is expressed by hormone-sensing LCs, triggering progesterone-dependent expansion of neighbouring Rank-positive cells, including BCs and LPs 4,9,10. Ubiquitous, whole mammary epithelial and LC-specific deletion of Rank impairs alveologenesis and lactogenesis due to compromised proliferation and differentiation of mammary epithelial cells 3-5. Gain-of-function mouse models, where Rank is mainly expressed by LCs (MMTV-Rank) show increased mammary epithelial cell proliferation and stemness, resulting in early alveologenesis and lactation failure 6-8. Thus, proper regulation of Rank signaling is crucial for MG epithelial development.

Besides, Rank signaling plays a significant role in mammary tumorigenesis: genetic or pharmacological inhibition of Rank signaling prevents or attenuates mammary tumorigenesis 7,11-17. Paradoxically, MMTV-driven Rank expression enhances LCs senescence delaying MG tumor onset, while promoting tumor aggressiveness once tumors are established 7.

Although Rank is expressed by BCs 4, BC-specific Rank signaling and its impact on MG development and cancer remain poorly understood. Our group has shown that luminal-specific deletion of Rank results in unfit LCs upon pregnancy, favoring the awakening of adult BC bipotency through Rank signaling, to replace incompetent LCs and restore lactation 3. Here, we demonstrate that, upon overactivation of Rank signaling, BCs reorganize their chromatin and undergo LC-like differentiation through epigenetic remodeling, leading to an aberrant basal-luminal hybrid cell population that eventually gives rise to mammary preinvasive lesions and ultimately breast adenocarcinomas. Moreover, we describe a basal Rank gene signature associated with recurrence of precancerous lesions and poor survival in breast adenocarcinomas.

Results
Enhanced basal Rank signaling disrupts basal and luminal cell identity
Our recent results suggest that Rank signaling plays a role in the postnatal reactivation of BC bipotency 3. To directly assess the role of Rank signaling in BCs, we generated a novel mouse model enabling the inducible expression of Rank specifically in BCs and allowing fate mapping of control (K14-Ctrl) or Rank (K14-Rank) overexpressing cells (Gfp+) at different developmental stages. Briefly, a doxycycline inducible Tet-On/Cre-loxP system was used to allow BC-specific Cre activation (under the control of keratin 14 (K14) promoter), subsequent Rank expression and tracing of recombined BCs and their progeny (see Methods for details) (Supplementary Fig. 1a). Experimental and control mice were treated with doxycycline during puberty (weeks 4-8) (Fig. 1a), and BC-specific recombination was confirmed by co-localization of Ck14 and Gfp. BC-Rank recombination remained stable over time in K14-Ctrl and K14-Rank MGs and co-localization of Gfp and overexpressed Rank in the BCs of K14-Rank MGs validated our model (Fig. 1b-d).

Whole mount analysis of MGs failed to demonstrate major differences in fat pad invasion or epithelial morphology upon BC-enhanced Rank expression compared to control females (Supplementary Fig. 1b). Immunohistopathological analysis determined a similar ratio of BC/LC and Hr+/Hr- LCs within the luminal epithelium (Fig. 1b, e and Supplementary Fig. 1c-f). BCs, but not LCs, initially showed enhanced proliferation in K14-Rank mice that was not sustained over time (Fig. 1f and Supplementary Fig. 1g, h). Unlike in the MMTV-Rank model 7, in which ectopic Rank expression is mainly targeted to LCs, no significant changes in the senescence marker SA-Gal were observed in BCs or LCs upon basal Rank expression (Supplementary Fig. 1i, j).

Next, we characterized the transcriptome of isolated Gfp+ BCs and Gfp- LCs from MGs of virgin K14-Ctrl and K14-Rank mice, and a complementary basal-specific Rank depletion mouse model (K14-RankΔ/Δ), in which no major alterations in MG development were detected 3. High levels of basal-specific recombination were confirmed in all genotypes (Supplementary Fig. 2a, b). LCs from all models clustered together upon PCA analysis, while BCs, particularly from K14-Rank mice, showed greater heterogeneity (Fig. 1g). Sample identity was confirmed by heatmap analysis of known luminal and basal markers (Supplementary Fig. 2c). K14-Rank BCs showed 211 differentially expressed genes (DEGs) compared to corresponding controls, including the upregulation of Mki67, in line with the increased levels of proliferation, but intriguingly, also LC-associated genes (Krt18, Ehf and Mfge8) and downregulation of myoepithelial/BC-associated genes (Oxtr, Fos and Egr1) (Fig. 1h, i and Supplementary Data 1).

No DEGs were found between BCs from K14-RankΔ/Δ and control littermates (Supplementary Fig. 2d), but gene set enrichment analysis (GSEA) revealed an opposite effect of Rank signaling in BCs according to Rank expression: K14-Rank BCs showed largely upregulation of gene sets (vs K14-Ctrl BCs), whereas in K14-RankΔ/Δ BCs most pathways were downregulated (Supplementary Fig. 2e and Supplementary Data 1). BCs from K14-Rank females upregulated gene sets associated with proliferation, Rank signaling (i.e. Nfb signaling and inflammation/immunity), protein translation and cell differentiation (i.e. Slits/Robo regulation), while most of these gene sets were downregulated in K14-RankΔ/Δ BCs, confirming Rank functionality in K14 driven gain and loss of function models, respectively (Supplementary Fig. 2f and Supplementary Data 1).

Considering the increased gene expression levels of luminal markers, we interrogated the K14-Rank BC gene data set against known BC and LC lineage differentiation gene signatures from the mouse and human breast 18; 19(Supplementary Data 2). In line with the role of basal Rank signaling modulating BC bipotency and luminal differentiation 3, BCs from K14-Rank mice show upregulation of signatures corresponding to mouse LP, mouse alveolar (AV), as well as human luminal AV cells, in particular the basal/luminal (BL) subpopulation (Fig. 1j and Supplementary Fig. 2g). Moreover, we computed a basal RANK signature (Rank-BC) considering the top upregulated (Rank-BC-UP) and downregulated (Rank-BC-DN) DEGs derived from the transcriptomic profile of K14-Rank BCs (vs K14-Ctrl BCs) (Supplementary Data 1, see Methods for details). Rank-BC gene signature correlated with luminal /alveolar-secretory progenitor (LASP) cells from healthy donors of the human breast cancer atlas (HBCA), including the LASP2 cluster characterized by mixed lineage expression 20 (Supplementary Fig. 2h). In these datasets, Rankl is expressed by hormone-sensing LCs (Supplementary Fig. 2i, j). Meanwhile K14-RankΔ/Δ BCs showed downregulation of gene signatures associated with both BC and LC lineages when compared to corresponding controls (Fig. 1k).

Then, we aimed to understand how the modulation of BC Rank signaling affects LCs. Gfp- LCs isolated from K14-Rank and K14-RankΔ/Δ showed none or few DEGs when compared to K14-Ctrl LCs, respectively (Supplementary Fig. 2k). Nonetheless, these Gfp- LCs shifted their identity when evaluated against known mammary epithelial lineage differentiation gene signatures 18 (Supplementary Fig. 2l).

These results demonstrate that the modulation of Rank signaling in BCs alters not only BC, but also LC identity, not reflected by an obvious phenotype in the MGs from virgin K14-RankΔ/Δ 3 and K14-Rank mice (Fig. 1b-f and Supplementary Fig. 1b-j).

BC-specific enhanced Rank signaling interferes with MG functionality
Considering the transcriptomic changes identified in BCs and LCs from virgin K14-Rank mice, we postulated that their MG functionality might be affected. To address the implications of ectopic basal Rank on lactation, we submitted K14-Rank females to a cycle of pregnancy and included untreated (UT) K14-Rank females to account for putative side effects related to K14-associated leakiness 3 (Fig. 2a). Near one third of the K14-Rank females lost their litter (Fig. 2b) and K14-Rank nursed-pups presented reduced weight compared to pups from K14-Ctrl and UT females, indicative of defective lactation (Fig. 2c). Whole-mount and histological analysis at L1 revealed less alveolar coverage in K14-Rank than in K14-Ctrl MGs, coupled with an expansion of BCs (Fig. 2d, e and Supplementary Fig. 3a, b). K14-Rank and K14-Ctrl MGs presented high levels of BC recombination at L1 (Fig. 2f, g). Intriguingly, L1 MGs from K14-Rank mice showed Gfp+ LCs, with higher Rank levels compared to K14-Ctrl mice (Fig. 2f, h, i and Supplementary Fig. 3c), in line with the gain of luminal gene expression markers by BCs and the shift to a luminal identity (Fig. 1i, j). Moreover, K14-Rank L1 MGs displayed an aberrant epithelial architecture characterized by small alveoli overexpressing Rank, similar to those found in MMTV-Rank females 8, and underdeveloped alveoli where only BCs overexpressed Rank (Fig. 2i). Further morphological characterization revealed the presence of preinvasive lesions exclusively in K14-Rank L1 MGs (Fig. 2j). These preneoplastic lesions were present in all K14-Rank L1 MGs and contained Rank+ cells invading the lumen, suggesting a BC origin (Fig. 2j, k). K14-Ctrl and UT K14-Rank mice presented similarly low percentages of Gfp+ BCs and LCs at L1 (Supplementary Fig. 3d), discarding a contribution of K14-related leakiness to the presence of Gfp+ LCs and observed phenotypes. These results demonstrate that enhanced BC-specific Rank signaling favors a discrete basal to luminal (B/L) transition during pregnancy and impairs LCs functionality leading to aberrant alveologenesis and compromised lactation.

Enhanced BC Rank signaling leads to BC-derived preinvasive lesions and the appearance of B/L hybrid cells
Pregnancy supports MG Rank signaling due to progesterone-driven Rankl expression. To confirm that enforced basal Rank signaling promotes the preinvasive lesions found at L1, we mimicked a pregnancy-like state inducing Rankl production by treating adult K14-Ctrl and K14-Rank females with medroxyprogesterone acetate (MPA), a synthetic progestin that acts as a potent agonist of progesterone receptor (Pr) (Fig. 3a). Analyses of MPA-treated K14-Ctrl and K14-Rank females confirmed Rankl overexpression and proto-alveolar development (Supplementary Fig. 3e, f), validating MPA activity. BCs from K14-Ctrl and K14-Rank MGs sustained full recombination upon MPA treatment (Supplementary Fig. 3g, h). MGs from MPA-treated K14-Rank females showed an expansion of the BC compartment coupled with enhanced BC, but not LC, proliferation (Supplementary Fig. 3i-l). Indeed, MPA treatment boosted the incidence of preinvasive lesions in K14-Rank MGs, with full penetrance as observed at L1 (Fig. 3b and Fig. 2j). Upon careful examination, similar preneoplastic lesions, reminiscent of ductal carcinoma in situ (DCIS), could also be detected in a few (2 out of 7) K14-Rank virgin mice (Supplementary Fig. 4a).

Immunohistochemistry (IHC) and immunofluorescence (IF) analyses revealed that these preneoplastic lesions originated from BCs, as they were positive for Rank and Gfp expression (Fig. 3b and Supplementary Fig. 4a-c). Cells within the lesions did no longer express BC/myoepithelial markers such as Ck5 and Sma1 (Supplementary Fig. 4d, e) but gained the expression of LC and LP markers (Ck8 and Ck6a) (Supplementary Fig. 4f), indicative of lineage identity dysregulation. Moreover, cells within the preinvasive lesions showed a B/L hybrid phenotype as they co-expressed Ck14, Ck8 and Gfp (Fig. 3c). Indeed, flow cytometry analysis of MPA-treated K14-Rank MGs unveiled a distinct Gfp+ B/L hybrid population (CD24hi, CD49fhi), absent in MPA-treated K14-Ctrl MGs (Fig. 3d). Altogether, our data indicate that, upon enhanced basal Rank signaling, BCs lose their BC-restricted lineage specification and generate B/L hybrid cells that could give rise to BC-derived intraductal preinvasive lesions.

Enhanced Rank signaling leads to BC-derived mammary tumor development
The loss of restricted cell identity has been recently associated with breast cancer risk 21,22. Since the preinvasive lesions formed by BC-derived Rank+ hybrid cells resembled DCIS, precursors of invasive breast cancer 23, we hypothesized that they would eventually progress to invasive adenocarcinomas. Indeed, enhanced basal Rank signaling led to spontaneous mammary tumors in aged (approx. 70-week-old) virgin K14-Rank females. Whereas all K14-Ctrl females died without palpable mammary tumors, 60% of K14-Rank females developed MG tumors (Fig. 3e, f). Aged K14-Rank females showed multiple BC-derived mammary intraductal preinvasive lesions in their MGs, harboring Gfp+ as well as Rank overexpressing cells and cells co-expressing basal Ck14 and luminal Ck8 markers (Supplementary Fig. 5a, b). Immunohistopathological characterization of aged K14-Rank spontaneous tumors indicated they were Er- adenocarcinomas, positive for Gfp and Rank overexpression, denoting their basal origin (Fig. 3g and Supplementary Fig. 5c). Like the aberrant B/L hybrid population in MPA treated mice (Fig. 3c and Supplementary Fig. 4d-f), K14-Rank BC-derived spontaneous tumors have lost the expression of some basal markers while gained luminal ones, being some tumor cells positive for both basal and luminal markers (Fig. 3h). Aiming to gain further insights into basal Rank induced tumorigenesis, we induced mammary carcinogenesis by using MPA and DMBA treatment 14,24 (Fig. 3i). K14-Rank females showed similar tumor incidence but accelerated mammary tumor onset compared to control mice (Fig. 3j, and Supplementary Fig. 6a). Histological characterization of MPA/DMBA-induced mammary tumors revealed the expected tumor type heterogeneity in the control group, with a 64.7% incidence of adenocarcinomas but also myoepitheliomas and squamous tumors, while K14-Rank females almost exclusively developed adenocarcinomas (91.3%) (Fig. 3k and Supplementary Fig. 6b). Half of the tumors from K14-Rank mice expressed estrogen receptor (Er) while K14-Ctrl tumors were predominantly Er+ (Supplementary Fig. 6c, d). Lineage tracing analysis revealed that most K14-Rank tumors (86.9%) derived from Gfp+ Rank overexpressing BCs; they were classified as adenocarcinomas (60% Er- and 40% Er+) and sustained high Rank expression levels (Fig. 3l-n and Supplementary Fig. 6e, f). Only 3 out of 23 tumors in K14-Rank females, all Er- adenocarcinomas, were Gfp-, indicative of a luminal origin (Supplementary Fig. 6g, h), and 1 out of 7 females developed a lung metastasis, not detected in any K14-Ctrl mice (Supplementary Fig. 6i, j). These results indicate that enhanced BC Rank signaling exacerbates BC oncogenic potential, accelerates mammary tumor development and aggressiveness, and favors the appearance of adenocarcinomas versus other tumor types.

Importantly, 46.6% of K14-Ctrl derived tumors (7 out of 15) were also Gfp+, therefore arising from BCs (Fig. 3l and Supplementary Fig. 6e). These Gfp+ tumors from control mice included Er+ adenocarcinomas and Er- adeno-squamous carcinomas (Fig. 3m, n), indicating that BCs with physiological levels of Rank generate luminal mammary tumors.

To assess whether endogenous basal Rank signaling is required for mammary tumorigenesis, we treated K14-RankΔ/Δ females with MPA/DMBA (Fig. 3i). The loss of basal Rank showed a protective effect towards mammary tumorigenesis as all K14-Ctrl mice developed tumors (4 out of 4), whereas only 60% of K14-RankΔ/Δ females did (3 out of 5) (Supplementary Fig. 7a). In addition, tumor latency was delayed in K14-RankΔ/Δ compared to control MGs (Supplementary Fig. 7b). Although some K14-RankΔ/Δ females developed adenocarcinomas (Supplementary Fig. 7c), they lacked Gfp+ tumor cells (Supplementary Fig. 7d), denoting that they were originated by LCs. These results demonstrate that BCs with physiological levels of Rank can instigate mammary tumorigenesis similar to LCs, and that loss of basal Rank signaling prevents mammary tumor onset. Besides, the loss of BC lineage fidelity induced by enhanced basal Rank signaling favors BC transformation and breast cancer development.

BC-derived hybrid cells show reorganization of chromatin accessibility consistent with aberrant luminal/alveolar differentiation
Mammary epithelial cell lineage commitment is controlled through the remodeling of the chromatin architecture. Repressive histone marks (such as H3K27me3 and methylated CpG islands) characterize inactive or closed chromatin, whereas active marks (H3K4me3, histone acetylation and unmethylated CpG islands) allow for open chromatin conformation and activation of gene transcription. Therefore, the epigenetic landscape of each cell lineage grants transcription of identity genes while turns off those characterizing other MG populations (reviewed in 25). To understand whether epigenetic perturbation caused the observed changes in BC/LC lineage identity, favoring the onset of a B/L hybrid population, we performed chromatin profiling experiments. ATAC-seq analyses of hybrid Gfp+ B/L cells, Gfp+ BCs and Gfp- LCs from MPA-treated K14-Rank and K14-Ctrl mice (Fig. 3d). ATAC-seq analyses revealed large differences in accessible chromatin between the K14-Rank B/L hybrid population and BCs/LCs from K14-Rank MGs (Fig. 4a, b and Supplementary Fig. 8a). The chromatin profile of known luminal and basal markers 26 in B/L hybrid cells from experimental females was in-between those of BCs and LCs, decreasing accessibility of BC-associated genes such as Krt5, Acta2 or Trp63 while increasing that of luminal ones (i.e. Krt6a, Krt8, Krt18 or Ehf) when compared to BCs and LCs (Fig. 4a, b, Supplementary Fig. 8a-c and Supplementary Data 3). Moreover, further examination of the ATAC-seq results unveiled a remarkable increased chromatin accessibility for the alveolar marker Csn2 in the B/L hybrid population compared to LCs and BCs (Fig. 4a, b and Supplementary Fig. 8d), suggesting that these intermediate cells also show an aberrant luminal/alveolar identity.

In regard to the representation of transcription factor (TF) motifs in the differentially accessible chromatin regions (considering unique ATAC-seq peaks), the B/L hybrid population showed substantial differences, with increased chromatin accessibility for binding sites of Fos/Jun, Batf3/Jun, and Nfb (Rel/a, Nfb2) TFs when compared to K14-Rank BCs and LCs, followed by enrichment for LC identity TFs Elf1-5 and Ehf (Fig. 4c, Supplementary Fig. 8e and Supplementary Data 3).

Of note, BCs and to a lesser extent LCs, also displayed differential chromatin accessibility profiles of known LC and BC-associated signatures 26. Rank-overexpressing BCs compared to BCs from controls show a clear upregulation of the luminal (LP/ML) chromatin profile (opening LP genes such as: Ehf, Trim46 and Wfdc18) (Fig. 4d, Supplementary Fig. 8f and Supplementary Data 3). Similar to hybrid B/L cells, motif enrichment for the TFs Batf3, Bnc2, Nfb (Rel/a) and Batf/Jun was detected in K14-Rank BCs relative to controls (Fig. 4e and Supplementary Data 3) in line with transcriptomic data (Supplementary Fig. 2f and Supplementary Data 1) and Rank-Nfb promoting BC bipotency 3. These results indicate that K14-Rank BCs remodeled their chromatin and acquired luminal features before undergoing a detectable phenotypic change and transitioning to the B/L hybrid population. Indeed, one of the biggest changes in gene chromatin accessibility when comparing K14-Rank BCs to control BCs corresponded to Ehf (Fig. 4d), whose mRNA was also found upregulated in the transcriptome of K14-Rank BCs (Fig. 1i). The TFs Ehf, RelA, Jun/Fos and Elf1-5 were recently described by Langille et al. 22 as facilitators of basal to luminal/alveolar-like lineage conversion and tumorigenesis upon loss of epigenetic regulators and histone modifiers such as Asxl2. To investigate how enhanced basal Rank signaling remodels epithelial protein programs, we performed quantitative TMT-based proteomic and phosphoproteomic analyses on freshly isolated epithelial organoids from virgin 8-week-old K14-Rank and K14-Ctrl females, which reflect the Rank gain-of-function epithelium. Rank-based proteomic profiling revealed that several chromatin regulators, including Kdm4a and Asxl2, together with CK5 occupied lower positions in the ranked distribution of K14-Rank samples relative to controls, whereas histone H1.3, Macro H2A1 and Ezh1, tended to occupy higher ranks (Fig. 4f, Supplementary Data 4). To interpret these shifts at the pathway level, we applied GSEA to the full ranked protein list and phosphosite rankings, which revealed enrichment of pathways associated with chromatin organization concomitant with a non-proliferative state (Supplementary Fig. 8g, Supplementary Data 4). These results suggest that activation of Rank signaling in BCs leads to chromatin remodeling and cell cycle arrest, which could be mediated through the regulation of epigenetic modifiers, resulting in basal to luminal conversion and generation of aberrant B/L hybrid cells.

Basal to luminal lineage conversion is altered by epigenomic targeting
Epigenetic reprogramming plays a critical role in the dynamics of the MG development and the differentiation of the mammary lineages (reviewed in 25). Given that our data pointed to chromatin remodeling as the mechanism underlying basal and luminal lineage infidelity in the context of basal Rank signaling, we then analyzed the global levels of the repressive mark H3K27me3 as it plays a central role in maintaining the mammary stem cells/progenitor pool and in mammary lineage commitment 27-30. K14-Rank BCs showed higher levels of H3K27me3 than K14-Ctrl BCs, in line with their luminal-like differentiation, whereas hybrid K14-Rank cells harbored increased H3K27me3 levels than BCs and lower than those detected in LCs (Supplementary Fig. 9a).

We then reasoned that Rank-driven regulation of H3K27me3 may contribute to B/L conversion and the appearance of the hybrid B/L K14-Rank cells. To assess this hypothesis, we made use of two epigenetic drugs, GSK-126 and GSK-J4, modifiers of H3K27 methylation status. GSK-126 suppresses H3K27me3 by inhibiting EZH2 methyltransferase, involved in maintaining luminal progenitors and restricting mature luminal differentiation 27,28. GSK-J4 is a H3K27me3 demethylase inhibitor targeting members of the KDM histone demethylase family, including KDM6 but also KDM4 31, resulting in the upregulation of H3K27me3 levels (Supplementary Fig. 9b). K14-Rank and K14-Ctrl females were treated with GSK-126 and GSK-J4 prior to stimulating basal Rank signaling with doxycycline and MPA (Fig. 4g). In the presence of GSK-126 or GSK-J4, K14-Rank and K14-Ctrl MGs showed decreased or increased H3K27me3 levels, respectively, compared with MGs from untreated mice (Supplementary Fig. 9c), confirming that both treatments were effective in vivo. Morphological assessment and flow cytometry analyses revealed no major alterations in K14-Ctrl MGs, with Gfp expression being exclusively detected in BCs (Supplementary Fig. 9d, e). However, in K14-Rank MG both drugs led to an increase in Gfp+ LCs (Fig. 4h, Supplementary Fig. 9e), suggesting that the modulation of H3K27 methylation levels promotes basal to luminal conversion. Importantly, the percentage of Gfp+ B/L hybrid cells was reduced in GSK-J4-treated K14-Rank MGs while the opposite trend was observed upon GSK-126 treatment (Fig. 4i). We next assessed the impact of these epigenetic inhibitors on the premalignant lesions found in K14-Rank MGs. In untreated K14-Rank mice, all preinvasive lesions contained Gfp+ B/L hybrid cells co-expressing Ck14 and Ck8 (Fig. 4j, k and Supplementary Fig. 9f). In contrast, a drastic reduction of lesions containing Ck14+Ck8+ hybrid cells was seen in GSK-J4-treated K14-Rank MGs where most lesions contained lineage-restricted Ck14+Ck8- BCs. GSK-126 treatment also led to a moderate decrease in the frequency of lesions containing hybrid cells (Fig. 4j, k). Altogether, these results confirm that disrupting the epigenetic landscape by altering H3K27me3 directly influences the B/L hybrid population endorsed by basal Rank signaling.

Basal Rank signaling is associated with DCIS recurrence and poor survival in luminal breast adenocarcinomas
Considering that the preinvasive lesions observed in K14-Rank MGs resemble DCIS, we investigated the clinical relevance of our findings in a collection of 57 DCIS and 313 invasive ductal carcinomas (IDC) 23. RANK and RANKL expression levels were comparable between DCIS and IDC samples (Supplementary Fig. 10a). When considering PAM50 classification, RANK expression was higher in basal IDC compared to other subtypes as previously reported 23,32,33, and a similar distribution of RANK expression was observed in DCIS, while RANKL was equally expressed between subtypes in DCIS and IDC (Supplementary Fig. 10b). The expression of the combined basal Rank signature (Rank-BC) was higher in IDC compared to DCIS (Fig. 5a) and higher levels of Rank-BC were found in basal, Her2 and LumB and the lowest in LumA IDC, whereas in DCIS the highest levels of Rank-BC were found in the Her2 subtype (Fig. 5b).

Importantly, a positive association between RANK and the transcription factor BATF3 was found in DCIS and IDC samples (Fig. 5c). Meanwhile, the association between RANK and ASXL2 or EZH2 followed opposite patterns in DCIS and IDC cases: in DCIS ASXL2 and EZH2 tended to be negatively associated with RANK whereas in IDC both genes were positively associated with RANK (Fig. 5c). Indeed, ASXL2 expression level was higher in DCIS than in IDC cases (Fig. 5d), which may indicate that the loss of function of the epigenetic modifier Asxl2, favored by Rank signaling, contributes to the progression from DCIS to IDC, in line with the preclinical findings. Besides, EZH2 and KDM6A showed opposite correlations with basalness, being EZH2 positively correlated and KDM6A negatively (Supplementary Fig. 10c), reinforcing their relevance in the modulation of basal/luminal identity in human breast adenocarcinomas.

Next, we tested whether our Rank-BC signatures predicted DCIS recurrence -an ipsilateral breast event (DCIS or IDC)- analyzing the independent DCIS panel Translational Breast Cancer Research Consortium (TBCRC) 038, which includes data of disease recurrence from 216 patients with clinically annotated longitudinal sample cohorts 34. Interestingly, the interrogation of the TBCRC 038 cohort revealed that Rank-BC-UP and Rank-BC-DN signatures correlated with DCIS recurrence (Fig. 5e and Supplementary Data 5). Moreover, pathways previously associated with DCIS recurrence (G2M/E2F pathways, MTORC1 signaling, and MYC targets V1) 34, were upregulated in K14-Rank derived BCs (Supplementary Data 1). Notably, the combined Rank-BC signature was predictive of any subsequent recurrence, DCIS or IDC (Fig. 5f), in line with the spontaneous progression of the precancerous lesions to invasive adenocarcinomas observed in the K14-Rank mouse model. Indeed, multivariate Cox regression analyses revealed that the basal Rank signature was predictive of recurrence (HR 2.02 (1.02-4.0) p=0.04), independent of age, clinical ER or DCIS grade (Fig. 5g).

Then, we evaluated the prognosis potential of the basal RANK predictive signature in IDC samples by interrogating the METABRIC dataset 35. Rank-BC signature was associated with poor survival in breast cancer patients from the METABRIC cohort, particularly in those harboring luminal breast tumors (Fig. 5h, i and Supplementary Fig. 10d). Together, our results denote that enhanced basal Rank signaling favors ipsilateral breast recurrence after DCIS, as well as poor survival in breast cancer patients, particularly in those with luminal breast cancer.

Methods
Generation of the knock-in B6.Cg-Gt(ROSA)26Sortm5(CAG-Rank,-GFP)Ics/N mouse line
The B6.Cg-Gt(ROSA)26Sortm5(CAG-Rank,-GFP)Ics/N mouse line was generated by inserting a floxed neomycin selection-STOP cassette followed by a construct containing murine Rank cDNA upstream of an IRES-eGfp-pA, in the Rosa26 locus. The targeting vector was electroporated into C57BL/6NCrl TB1 mouse embryonic stem cells, and clones were summited to neomycin selection (G418). Homologous recombinant clones were identified by Southern blot analysis of the genomic DNA digested with the restriction enzymes EcoRV and HindIII (neomycin probe I), Pacl and Xbal (neomycin probe II), and Xcml and SexAI (5’ probe); and using an external probe amplified from the neomycin cassette or from the genomic DNA. The band sizes obtained for the KI band for the neomycin probe I were 10.4 kb (EcoRV) and 8.0 kb (HindIII), for the neomycin probe II, 10.8 kb (Pacl) and 4.6 kb (Xbal) and for the 5’ probe, Xcml (14.1 (KI) and 6.0 (WT)) and SexAI (15.0 (KI) and 6.8 (WT)). Chimeras were generated by microinjection into BALB/cN embryos. Germline male chimeras were crossed with C57BL/6NCrl females to investigate whether the recombined ES cells have contributed to the germ layer. B6.Cg-Gt(ROSA)26Sortm5(CAG-Rank,-GFP)Ics/N were fertile and did not exhibit developmental phenotypes.

Mice generation and in vivo treatment
Animal husbandry, ethical handling of mice, and all animal work were carried at Spanish National Cancer Research Center (CNIO) and IDIBELL animal facility according to guidelines approved by the Committee on Animal Care and following Spanish and European Union regulations. All animal experiments were approved by the ethical Committee for Animal Experimentation (#CBA07_2021) from the Autonomous Community of Madrid (CAM). All mice used in this study were maintained in pathogen free cages under controlled condition of temperature (20-26ºC), humidity (30-70%), 12-hour light/dark cycles, food and water ab libitum. Chow diet was used 2018 Teklad Global 18% Protein Rodent Diets (Inotiv; https://www.inotiv.com/rodent-natural-ingredient-2018-diets). All mice used were females in C57BL/6JOlaHsd or mixed background. Control (K14-Ctrl) mice were: Tg.KRT14-rtTA: Tg.tetO-CMV-SCre (provided by Dr. Cedric Blanpain) and Rosa26_ACTB-tdTomato_EGFP (mRTmG, MGI: 3716464). Experimental mice (K14-Rank) carried previous allele combined with the R26Sortm5(CAG-Rank,-GFP)Ics/N allele were generated by C Mueller. K14-RankΔ/Δ mice were described in 3. Control and experimental mice were treated with doxycycline (1 mg/ml doxycycline (ITW Reagents, #A2951.0025) in 20 mg/ml sucrose (Merk, # S9378)), ad libitum in the drinking water to achieve Cre recombination, as described in 3. For enhancing Rankl expression, 6- or 8-week-old mice were subcutaneously treated with 50 mg/pellet of methyl-progesterone (MPA) (Innovative Research of America #NP-161) for 21 days. For modulating the epigenome (by targeting H3K27 methylation), 4-week-old mice were treated intraperitoneally 3 times per week during 5 weeks with 50 mg/kg of the EZH2 inhibitor GSK-126 (MedChemExpress, HY-13470-50 MG) or the KDM6 inhibitor GSK-J4 (MedChemExpress, HY-15648B-50 MG) dissolved in DMSO/Captisol® (Selleck Chemicals, S4592-5G). Mammary tumors were induced by MPA and 7,12-dimethylbenzanthracene (DMBA) (Fisher Scientific #408181000) as described in 14,24. Mice were monitored for tumor onset by manual palpation biweekly and measured by caliper once palpable tumors (3mm diameter) were detected. Mice were sacrificed once tumors reached 10 mm diameter or at humane endpoint. Maximum tumor size was not exceeded in this study.

Whole mount staining, immunofluorescence (IF) and immunohistochemistry (IHC)
MGs were excised and stretched onto a glass slide. For whole mount staining MGs were Carnoy fixed (ethanol 95%: chloroform: glacial acetic acid (6:3:1)) for 2h at room temperature (RT). MGs were washed with 70% EtOH and distilled water and incubated in carmine aluminum solution O/N at 4ºC. MGs were dehydrated by increased EtOH washes and stored in xylene. Visualization and capture were performed using a Leica MZ10-F stereo microscope. For IF and IHC, mouse tissue samples were formalin-fixed and embedded in paraffin. For IF, antigen heat retrieval (citrate buffer pH6) was used on paraffin-embedded sections unless indicated. 5 µm tissue sections were blocked with 10% BSA-5% donkey/goat serum prior primary antibody incubation O/N at 4ºC (primary antibodies: Ck8 (Troma-I, DSHB, 1:400, clone Troma-I), Ck14 (AF-64, Covance, 1:400, polyclonal), Ki-67 (556003, Biolegend, 1:300, clone B56-RUO), Ck5 (905901, Biolegend, 1:500, polyclonal), Sma-1 (A2547-100UL, Merk, 1:400, clone 1A4), progesterone receptor (Pr) (12683667, Fisher Scientific, manufacturer’s recommendations, clone SP2), Gfp (ab13970, Abcam, 1:1000, polyclonal), and Ck6a (PRB-169P, Covance, 1:50, polyclonal). Tissue sections were then incubated with corresponding fluorochrome-conjugated secondary antibodies: Alexa Fluor 488 Donkey anti rabbit (A-21206, ThermoFisher Scientific, 1:400, polyclonal), Alexa Fluor 555 Goat anti-rat (A21434, ThermoFisher Scientific,1:400, polyclonal), Alexa Fluor 647 Goat anti-rabbit (A21245, ThermoFisher Scientific, 1:400, polyclonal), Alexa Fluor 647-conjugated AffiniPure Donkey anti-chicken IgY (H+L) (A78948, ThermoFisher Scientific, 1:400, polyclonal), and Alexa Fluor 488 Goat anti-mouse (A11001, ThermoFisher Scientific, 1:400, polyclonal) and Dapi (Merk Life, #D9542, 1:400) to stain cell nuclei for 2 hours at RT and mounted with Prolong Gold Antifade reagent (Life Technologies #P36930). Confocal analysis was carried out using a Leica TCS SP5-MP confocal microscope. For IHC, tissue sections were incubated O/N at 4°C with primary antibodies (Rankl (462, RnD, 1:200, polyclonal), Rank (AF692, RnD, 1:200, polyclonal), Gfp (ab13970, Abcam, 1:1000, polyclonal) and Estrogen receptor (Er) (AM (CRET94D/H9), CNIO Monoclonal Antibodies Core Unit, manufacturer's recommendations, clone CRET94D/H9)) and detected with corresponding biotinylated-streptavidin horseradish peroxidase conjugated secondary antibodies from VECTASTAIN ELITE ABC kit Peroxidase (Vector laboratories, Manufacturer’s recommendations) and revealed with DAB substrate (Vector laboratories, #SK-4105). Sections were counterstained with hematoxylin. For Rankl IHC, tissue sections were treated with 100 μg/mL Pronase (Sigma-Aldrich, #10165921001) following antigen heat retrieval (citrate buffer pH6). For Rank, antigen retrieval was performed by tissue sections incubation with Protease XXIV (Sigma-Aldrich #P8038-100MG). For IF analysis, B/L hybrid cells are considered as Ck8+Ck14+ (co-expression of both markers), while luminal cells are identified based on Ck8+ and basal cells as Ck14+/Ck5+.

Epithelial cell isolation and flow cytometry sorting
Fresh MGs were mechanically dissected by scissor chopping after draining lymph nodes were removed and enzymatically digested with digestium medium (Dulbecco’s modified Eagle’s medium (DMEM)/F-12 (Biowest, #L0093-500ml), 0.6 U/ml collagenase A (Sigma-Aldrich, #1108879300), 2.85 U/mL dispase II (ThermoFisher Scientific, #17105041), 2% HEPES (ThermoFisher Scientific, #H8357-100 ml) and 1% penicillin–streptomycin (ThermoFisher Scientific, #15140-122)) for 1h rocking (50 rpm) at 37 °C. Samples were washed with CO2 independent media (ThermoFisher Scientific, # 18045088) supplemented with 2mM glutamax (Gibco, #35050-061). Single cells were isolated by sequential treatments with 0.25%trypsin/0.1% EDTA solution (Sigma-Aldrich, #T3924-100ml) for 1 min at RT and dispaseII/DNAseI (ThermoFisher Scientific, #17105041 / Roche, #10104159001) at 37°C for 5 min. Enzymatic digestions were blocked by adding CO2 independent media (ThermoFisher Scientific, #18045088) supplemented with 2mM glutamax (Gibco, #35050-061) and 5% FBS (ThermoFisher Scientific, #F7524-500ml). Erythrocytes were eliminated by treating samples with cold (4°C) ammonium chloride solution (STEMCELL Technologies, #7008) for 1min at RT. Cell aggregates were removed by filtering the cell suspension with a 40 μm strainer and single epithelial cells were counted manually using a Neubauer chamber. Single cells were blocked with blocking buffer (PBS + 0.5% Bovine Serum Albumin (Sigma-Aldrich, #13391+ 2mM EDTA) for 20 min on ice. Single cells were incubated for 40 min on ice with the corresponding surface antibodies (CD24-AlexaFluor700 (564237, BD Pharmingen™, 1 g/ml, clone M1/69), CD49f-AlexaFluor647 (562473, BD Pharmingen™, 1.25 g/ml, clone GoH3), CD31-PECy7 (102417, Biolegend, 0.25 µg/ml, clone 390), and CD45-PECy7 (103113, Biolegend, 0.125 µg/ml, clone 30-F11)). Prior sorting, single cells solution was filtered by a 70 μm strainer. Cell sorting was performed using a BD FACSAriaTM cell sorter. The markers to specify BC, LC and B/L hybrid cells are based on the surface expression of CD24 and CD49f, commonly used for flow cytometry characterization and FACS isolation. BCs are identified as: CD24Low, CD49fHigh; LCs as: CD24High, CD49fLow; and B/L hybrid cells as: CD24high, CD49High.

Bulk-RNA sequencing
Sorted cells were directly collected into the lysis buffer from RNeasy micro kit (Qiagen, #74004). Samples were mechanically homogenized by passing through 30G syringe 5-times, and RNA extracted according to the manufacturer’s instructions. From total RNA samples, 500 ng with an average RIN= 7.4 (range: 6.1–8.4) were processed into cDNA sequencing libraries with the “QuantSeq 3’mRNA-Seq Library Prep Kit (FWD) for Illumina” (Lexogen, #015) with a UMI Second Strand Synthesis module. Briefly, library generation is initiated by reverse transcription with oligodT priming, followed by a random-primed second strand synthesis. While primers from both steps contain Illumina-compatible sequences, random primers additionally feature 6 nt long Unique Molecular Identifier (UMI) tags. Libraries are completed by PCR and sequenced on an Illumina NextSeq 550 (with v2.5 reagent kits) by following the manufacturer’s protocols. This library preparation kit generates directional libraries stranded in the sense orientation, i.e., the read1 (the only read in single read format) has the sense orientation. RNA-seq data has been submitted to GEO ID: GSE272902.

Data processing and normalization
For the analysis of the RNA-seq reads, Nextpresso v2 55 was used, as follows: UMI sequences within reads were extracted with UMItools 56. The resulting FASTQ 57 files were processed with bbduk (https://sourceforge.net/projects/bbmap/) for the extraction of adapters, polyA and rRNA sequences. Read quality was verified with FastQC v0.11.7 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The derived FASTQ files were aligned to the mouse genome (GRCm39) with Hisat2 58. The obtained BAM files were sorted with SAMtools 59, and the UMItools were employed once again for read deduplication. Htseq-count 60 was used to count the number of reads that aligned within genes, using the Gencode vM29 gene annotation 61. ComBat-seq 62 was then used to adjust the batch effect.

Gene expression analysis
Differential expression was performed with DESeq2 63, keeping only those genes with more than 10 counts across the samples. Finally, those genes that had an adjusted p-value below 0.05 were selected. GSEAStandard 64 was used to perform gene set enrichment analysis (GSEA) for the selected gene signatures, setting 1000 gene set permutations. Only those gene sets with significant enrichment levels (FDR < 0.25) were considered. A minimal and maximum number of genes per reference were defined as 10 and 500 respectively. GSEA was performed with gene signatures identified by previous literature at a single-cell level 18. Summary GSEA bubble heatmaps were generated with ggseabubble 65.

ATAC-seq sample preparation and sequencing
For ATAC-seq, 25.000 FACs-isolated cells were washed in cold PBS (500g, 5’, 4ºC) and lysed in 50ul lysis buffer (10 mM TrisHCl pH7,5, 10 mM NaCl, 3 mM MgCl2, 0,1% NP40, 0,1% Tween 20, 0,01% Digitonin and 1x Complete Protease Inhibitor Cocktail) for 5 min on ice. Nuclei were isolated by centrifugation (500g, 10 min, 4ºC) and treated with transposase Tn5 (Illumina, Nextera DNA Library Preparation Kit, #FC-131-1024). DNA was purified using MinElute PCR Purification Kit (Qiagen, #28004). Optimal cycles for library preparation were determined by agarose gel based on nucleosome pattern integrity. All samples were then amplified by PCR using NEBNextHigh-Fidelity 2× PCR Master Mix (New Englands Labs, #M0541S) with primers containing a barcode to generate libraries. Subsequently, size selection of DNA fragments was performed using Sera-Mag Select Beads (Thermo Fisher, #09-981-R123). A 0.7x volume of beads mixture was used to remove fragments larger than 600 bp, followed by the addition of an additional 1.1x volume to isolate fragments larger than 100 bp. Purified libraries were sequenced using NovaSeq system (Illumina) Pair End Run (64 and 150 bp) type at the Genomic Unit of CNIO (Madrid, Spain). ATAC-seq data has been submitted to GEO ID: GSE272968.

ATAC-seq data processing
Sequenced reads obtained from NGS were demultiplexed and quality filtered using FastQC (Andrews, 2010). Subsequently, Bowtie 2 66 was used for paired-end read alignment against mouse (mm9) genome. Signal tracks were generated with RPKM normalization, using a bin size of 50 bp. For quantification of ATAC-seq data, differential signal analysis was performed at ENSEMBL annotated genes through the integration of R libraries, including 'csaw', 'edgeR', 'rtracklayer', and 'statmod'. Bam files were used as input with the specific parameters: window width = 150, spacing = 50 and minq = 40. ‘normFactors’, ‘asDGEList’, and ‘estimateDisp’ functions were used for normalization and differential binding analysis using CPM values. Subsequently, the ‘glmQLFit’ and ‘glmLRT’ functions were used for comparative tests.

Visualization of signal tracks was performed using the UCSC Genome Browser, with mean values and a smoothing window between 6-12 pixels. MACS2 67 was used for calling broad peaks from ATAC-seq experiments. Peaks were result of common peaks between replicates. Differentially enriched motifs were identified using the Sequence Enrichment Analysis (SEA) tool based on the unique peaks in each comparison. The SEA tool (v5.5.5) was run using default parameters comparing sequences from ATAC-seq peaks against control sequences (shuffled input sequences). Known motifs from the JASPAR CORE 2022 vertebrate database were scanned. Enrich TF criteria (enrichment > 1, pval ≤ 0.05).

Epithelial organoid harvesting
Fresh MGs were mechanically dissected by scissor chopping after draining lymph nodes were removed and enzymatically digested with digestium solution (Dulbecco’s modified Eagle’s medium (DMEM)/F-12 (Biowest, #L0093-500ml), 2 mg/ml collagenase A (Sigma-Aldrich #1108879300), 2 mg/mL Trypsin (ThermoFisher, #27250018), 5% FBS (ThermoFisher Scientific, #F7524-500m), 5mg/ml insulin (Sigma-Aldrich, #I6634-50MG) and 50 mg/ml gentamicin (Sigma-Aldrich, #G1272-10MG) for 30 min rocking (55 rpm) at 37 °C. Samples were washed twice with DMEM/F-12 (Biowest, #L0093-500ml). Organoids were isolated by sequential treatments with Red Blood Cell Lysing Buffer Hybri-Max (Sigma-Aldrich, #R7757-100ML), for 1 min at RT, 20U/ml DNAseI/ 0.5mg/ml dispase II (Roche, #10104159001 / ThermoFisher Scientific, #17105041) at 37°C for 5 min, and 5 rounds of differential centrifugations (450g for 10 sec). Pelleted organoids were collected in Protein LoBind 1.5ml tubes (Eppendorf, #022431081) and snap-frozen at -80ºC before lysis for proteomic analysis.

Sample preparation for proteomic analysis
As FACS-sorted basal or luminal cells did not yield sufficient protein for deep mass spectrometry analysis, unbiased proteomic profiling was performed on freshly isolated organoids, which constitute an enriched epithelia fraction. Epithelial organoids were isolated from virgin 8-week-old females and processed immediately after isolation. For this study, five independent K14-Ctrl and five independent K14-Rank biological replicates were analyzed. As no in vitro stimulation, culture, or treatment was applied at any point, the resulting proteomic and phosphoproteomic datasets reflect the basal, untreated in vivo epithelial state under Rank gain-of-function conditions. Following isolation, organoids were directly lysed for 15 minutes at 95°C in 4% SDS, 100 mM Tris-HCl (pH 8.0). After cooling, lysates were treated with 10 units of DNase (Benzonase, Merck) at 25°C and sonicated for 10 minutes using a Bioruptor to shear DNA. Protein concentration was quantified using a microBCA assay with BSA as the standard. Protein digestion was performed using the on-bead Protein Aggregation Capture (PAC) method with MagReSyn® Hydroxyl microparticles (protein-to-beads ratio 1:5) on a King Fisher automated instrument (ThermoFisher). Proteins were digested overnight (16 hours) at 37°C with 300 µL of trypsin/LysC in 50 mM TEAB (pH 8.0) at a protein-to-enzyme ratio of 1:100 (Trypzean, Sigma; LysC, Wako). The resulting peptides were dried in a speed vacuum and re-dissolved in 100 µL of 200 mM HEPES (pH 8.5).

TMTpro Labeling
Samples (~100 µg) were labeled for 1 hour at 25°C using Thermo Scientific TMTpro 18plex™ isobaric reagents following the manufacturer’s instructions. Each of the 18 TMT channels corresponded to an independent biological replicate, derived from distinct mice. TMTpro labeling efficiency was evaluated using a pooled sample. After labeling, reactions were quenched with 5% hydroxylamine. Labeled samples were combined in a 1:1 ratio based on total peptide amounts, determined by comparing overall signal intensities in an LC-MS/MS run. The pooled mixture was desalted using a Sep-Pak C18 cartridge (Waters) and dried before high-pH reverse-phase HPLC pre-fractionation. Although a full TMTpro 18-plex design was used to analyze 18 fully isogenic samples, only the data corresponding to TMT channels 126–131C are reported in this study. For confidentiality reasons, the authors have decided to restrict the scope of this manuscript to this analysis and display of these samples exclusively. The remaining samples, although part of the same TMT labeling and MS acquisition run, are not presented here and will be described separately.

High-pH Reverse-Phase Chromatography
Peptides were pre-fractionated offline using an Ultimate 3000 HPLC system equipped with a sample collector. Peptides were dissolved in 100 µL of phase A (10 mM NH₄OH) and loaded onto an XBridge BEH130 C18 column (3.5 µm, 150 mm length, 1 mm ID, Waters). Phase B consisted of 10 mM NH₄OH in 90% acetonitrile (CH₃CN). Separation was performed using a gradient at a flow rate of 100 µL/min:0–50 min: 0–25% B, 50–56 min: 25–60% B and 56–57 min: 60–90% B. One-minute fractions were collected from 15 to 65 minutes, neutralized with 10 µL of 10% formic acid, and vacuum-dried. Based on UV absorbance at 280 nm, 40 fractions were combined into 8 pools for phosphopeptide enrichment.

Phosphopeptide Enrichment
Phosphopeptides were enriched using MagReSyn® Zr-IMAC HP beads on a King Fisher automated instrument following the manufacturer’s protocol. Eluted fractions were acidified with 10% formic acid and vacuum dried. Flow-throughs from each pool underwent further fractionation into four subfractions using micro RP-high-pH chromatography and were reserved for total proteome analysis.

Mass Spectrometry
LC-MS/MS analysis was performed on an Orbitrap Exploris 480 mass spectrometer coupled to an UltiMate 3000 RSLCnano system. Samples were loaded onto a trap column (Acclaim™ PepMap™ 100 C18, 5 µm, 20 mm length) at 10 µL/min in 0.1% formic acid (FA) for 3 minutes. Peptides were then separated on an EASY-Spray PepMap RSLC C18 column (2 µm, 75 µm × 50 cm, Thermo) at 45°C using a 90-minute gradient at 250 nL/min:0–2 min: 4–6% buffer B (100% ACN, 0.1% FA), 2–60 min: 6–33% buffer B and 60–70 min: 33–98% buffer B. The mass spectrometer was operated in a data-dependent mode with a top-15 acquisition method. MS spectra were acquired at 60,000 resolution (200 m/z), covering 350–1500 m/z. Peptide ions were isolated with a 0.7 Th window and fragmented via higher-energy collisional dissociation (HCD) at a normalized collision energy (NCE) of 36. MS/MS spectra were acquired with a fixed first mass of 120 m/z at 45,000 resolution. Ion targets were 3e6 for MS scans (maximum IT: 25 ms) and 1e5 for MS/MS scans (maximum IT: auto). The mass spectrometry proteomics and phosphoproteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE partner repository with the dataset identifier PXD054102 (proteomics) and PXD058134 (phosphoproteomics).

Proteomic and Phosphoproteomic Data Analysis
Raw files were processed using MaxQuant (v2.1.4.0) with TMTpro 18-plex quantification enabled. Database searches were performed against the UniProtKB/TrEMBL mouse reference proteome (21,990 sequences; downloaded March 1, 2023) supplemented with a curated contaminant list routinely used in the CNIO Proteomics Unit. Carbamidomethylation of cysteine and TMTpro (+304.2071 Da) on peptide N-termini and lysine residues were set as fixed modifications. Oxidation of methionine, N-terminal acetylation and deamidation (N/Q) were set as variable modifications. A minimum peptide length of seven amino acids and up to two missed cleavages were allowed. Peptide-spectrum matches, proteins and phosphosites were filtered at a 1% false discovery rate (FDR). TMTpro labelling efficiency and potential over-labelling were assessed via dedicated MaxQuant and Proteome Discoverer searches of pooled TMT samples. Labelling efficiency exceeded 98%, with <1% of total precursor intensity corresponding to unlabelled peptides. Over-labelling at serine, threonine and tyrosine residues accounted for 4.99% of total precursor intensity (6.46% of PSMs), well below reported values for TMTpro workflows, indicating minimal impact on quantitative accuracy. Reporter ion intensities were log₂-transformed and processed using Prostar (v1.30.0) and DAPAR (v1.22.3). LOESS normalization was applied independently to the proteome and phosphoproteome datasets. Entries with fewer than 18 valid reporter values were excluded. Together, these metrics confirm that TMTpro labelling was highly efficient with minimal unintended STY modifications, supporting the robustness of our quantitative proteomics workflow.

Because TMT-based quantification at the MS2 level (without FAIMS) is prone to reporter ion interference and ratio compression, and therefore not optimized for estimating small effect sizes at the individual-protein level 68, we did not apply fixed fold-change cutoffs or multiple-testing–adjusted significance thresholds to define sets of differentially expressed proteins. Although per-protein statistical scores (e.g., p-values) are reported for visualization (e.g., volcano plots), downstream analyses were based exclusively on rank-based approaches that capture coordinated pathway-level trends. For the proteome, proteins were ranked solely by log₂ fold change (K14-Rank vs K14-Ctrl) and used as input for gene set enrichment analysis (GSEA; v4.2.3) with the mouse msigdb_v2022.1.Mm collection (GMT format). For phosphoproteome-level analysis, GSEA was performed independently of protein abundance changes. Quantified phosphosites were first collapsed to a single value per gene to enable gene-centric enrichment analysis. For each gene, the phosphosite exhibiting the maximum absolute log₂ fold change across conditions was selected, and the signed log₂ fold change of that site was used to rank genes. This strategy captures the dominant phosphorylation change per protein while avoiding dilution or cancellation effects arising from phosphosites with opposing regulation on the same protein.

Protein-level fold changes were not used to normalize, weight, or adjust phosphoproteomic rankings, even when the corresponding protein was detected in the global proteome dataset. Proteome and phosphoproteome GSEA analyses were therefore conducted as independent analyses reflecting changes in protein abundance and phosphorylation state, respectively. All GSEA analyses were performed using the full ranked lists and the complete set of quantified proteins or phosphoproteins as background, consistent with a rank-based framework designed to detect coordinated pathway-level shifts rather than individual differentially expressed entities

Generation of BC and LC predictive Rank associated gene signatures
Rank associated gene signatures from BCs was generated from transcriptomic data of K14-Rank BCs (vs K14-Ctrl BCs) (Supplementary Data 1). Genes were clustered in upregulated and downregulated (FC>0; FC<0). Clustered genes were sorted based on adjusted p-value. The top 100 upregulated and downregulated genes with smallest p-value were considered associated Rank-BC-UP and Rank-BC-DN gene signature. Basal RANK signature (Rank-BC) was generated combining Rank-BC-UP and Rank-BC-DN gene signatures. The BC-DN genes were included as negative.

Public data analysis
For the analyses of the Human Breast Atlas 20 the tl.score_genes function in Scanpy calculates a score for a specified set of genes by computing the average expression of these genes and subtracting the average expression of a reference set. This reference set is randomly sampled from the gene pool, ensuring that the score reflects the relative expression of the target genes compared to a background 69.

The METABRIC dataset 35 (n = 2,509) was downloaded from cBioPortal (http://www.cbioportal.org/), and TBCRC data 34 was obtained from Human Tumor Atlas Network (HTAN, https://humantumoratlas.org). K14-Rank BC signatures (up- or down-regulated DEGs vs. K14-Ctrl BCs) were used to assign a score to each patient using single-sample GSEA (ssGSEA) 70. BC-DN score was subtracted from BC-UP score, and a global BC signature score was used. Patients were then ranked based on their signature score, and highest and lowest quartiles (Q1 and Q3) were compared for the clinical outcome of interest (survival or recurrence). Kaplan-Meier survival curves were used to visualize the differences in outcome, and association of Rank signature with time to recurrence or death was quantified by Cox Proportional Hazard model.

Gene expression analyses of DCIS and IDC
Previously published gene expression data from 57 DCIS and 313 IDC were obtained and PAM50 subtyping was performed as described in 23. Signature Z-scores were calculated as the mean of the standardized gene expression of all genes within the BC-UP and BC-DN signatures present in the dataset. The BC-DN genes were included as negative numbers. Correlation between gene centred expression data of the PAM50 genes and the basal-like centroid was calculated for each tumor sample using Spearman correlation as described in 23,33.

Statistics & reproducibility
Statistical analyses were performed using GraphPad Prism software version 10. Data are represented as the mean +/- SEM. The measurements were taken from different mice, and we avoided technical replicates. Statistical analyses are described in figure legends. All investigators were blind to the time-point and treatment arms. Quantification of IF/IHC and data analysis was performed without knowledge of the identity of the samples.

Data Availability
The raw data generated for RNA-seq analysis in this study have been deposited in the Gene Expression Omnibus (GEO) database under the ID GSE272902 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE272902]. The mouse genome GRCm39 was used for the transcriptomic analyses. The ATAQ-seq raw data has been deposited in GEO under the ID: GSE272968 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE272968] and alignment to the genome can be found here [http://genome.ucsc.edu/s/JoseTerron/mm9_jaime]. The proteomic and phospho-proteomics data have been deposited to the ProteomeXchange Consortium with the dataset identifier PXD054102 [https://www.ebi.ac.uk/pride/archive/projects/PXD054102] and PXD058134 [https://www.ebi.ac.uk/pride/archive/projects/PXD058134]. Data is publicly available. Biological materials can be obtained upon request. Source Data are provided with this paper.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript describes an investigation of the role of Rank signaling in basal cells of the mammary gland. The authors developed a novel mouse model for the inducible expression of Rank specifically in basal cells which allows for the lineage tracing of Rank-overexpressing cells at different developmental stages. Here, they showed that enhanced Rank signaling in basal cells resulted in the reorganization of their chromatin accessibility, favoring basal-luminal conversion which leads to mammary preinvasive lesions and ultimately breast adenocarcinomas.\nOverall, I found the manuscript well written, with clear description of the methods and data figures, and their observations of the role of Rank signaling in basal cell identity and tumorigenesis very intriguing. However, I have some key concerns regarding how this signaling pathway is activated in breast cancer patients, the major regulators underlying these chromatin accessibility changes and the key effector genes leading to the basal-luminal cell conversion and ultimately breast carcinogenesis. It will be helpful if the authors can clarify these to enhance the clinical significance of their studies, mainly conducted in mouse models, and address the following comments in their revised manuscript.\n1. Although the authors showed that simulating pregnancy via MPA can induce Rankl overexpression, they should clarify what are the natural sources or cell types producing Rankl in the mammary gland and corroborate this with expression data from human DCIS or luminal breast cancer patients, either using scRNA-seq or IHC analysis.\n2. The drastic changes in chromatin accessibility upon enhanced Rank signaling is a very interesting finding. While the authors demonstrated these chromatin remodeling changes are associated with various transcription factors and epigenetic regulators like Asxl2, they should verify which of these are the critical regulator(s) of this chromatin reorganization. This can be validated in vitro using mouse/human mammary or breast cancer cell lines through targeting these transcription factors and epigenetic regulators under enhanced Rank signaling.\n3. Additionally, the clinical relevance of their findings can be enhanced with further confirmation of whether these chromatin remodeling changes and Rank transcriptional signature can be similarly induced in human invasive or luminal breast cancer cells upon Rankl stimulation.\n4. Finally, the authors should perform functional studies of the target genes of the key regulator(s) found to underscore these chromatin remodeling changes to explain their biological roles in mediating the basal-luminal cell conversion and tumorigenesis.", "reviewer_2": "The study by Redondo-Pedraza et al. used a lineage-specific expression of the Rank model (K14-driven for specific BC expression) with dox-inducible expression at a specific time of development (puberty) to demonstrate that basal cell-specific RANK expression disrupts the distinct identities of basal and luminal mammary epithelial cells. Chronic activation of RANK in basal cells diminishes key basal lineage markers while aberrantly inducing luminal differentiation markers, resulting in lineage bipotency. They show that this defect in lineage commitment impairs normal mammary gland function, including lactation, and leads to the emergence of premalignant lesions composed of basal-derived cells with mixed luminal and alveolar features, which can progress to both basal and luminal breast adenocarcinomas. They that RANK activation could drive the loss of tumor-suppressive epigenetic regulators (AsxL2, Kdm4a), leading to chromatin remodeling, increased access to the repressor Batf3, and facilitating tumorigenesis. The study also identified a basal cell RANK gene signature associated with poor prognosis (in luminal adenocarcinomas) and recurrence in breast cancer, suggesting the clinical significance of basal RANK signaling in early breast cancer development and progression.\nOverall, the study uses an interesting and robust system for BC-specific and inducible RANK overexpression to demonstrate that aberrant RANK expression in BC impairs lineage specificity and contributes to lactation defects and pre-tumorigenic lesions, contributing to breast adenocarcinoma. These observations are interesting and could be of significance, especially because BC cells express RANK signaling. However, some conclusions are not fully supported by the data. The study lacks several important controls, and some key aspects remain to be addressed to mechanistically show that RANK signaling contributes to these effects and that this is occurring through the deregulation of epigenetic modulators and profiles.\nMajor comments:\n1- The study states that “RANK signaling” is responsible for basal cell-lineage infidelity, leading to mammary gland dysfunction, problems in lactation, and tumorigenesis. However, there is no confirmation of the functional activation of RANK signaling in the K14-driven RANK overexpression model compared to the Ctl model (e.g., downstream activation of RANK effectors). While the overexpression of RANK should likely lead to increased downstream RANK signaling, this should be shown (and perturbed) to ensure that the observed phenotypes and conclusions are directly caused by the activation of RANK signaling in the model (especially that BC cells already display some levels of RANK signaling, see point #2 below).\n2- It is stated in the introduction that BC expresses RANK signaling. This baseline level of RANK signaling in BC should be shown (again via showing activation of downstream RANK signaling effectors) and compared to the dox-induced overexpression of RANK in BC cells. Is the RANK signal induced by progesterone or MPA in the K14-Clt cells? And how does it compare to the aberrant signaling of the K14-RANK cells? Importantly, whether the aberrant RANK signaling resulting from RANK overexpression is relevant physiologically should be addressed.\n3- Some conclusions are based on tendencies (i.e. Not significant changes but with a trend, e.g. S3a, 1f, etc.). However, in other cases, this type of trend is claimed non-significant but is mainly due to large variability of the data points (eg. 1e,f,, 2e, S1j, etc). p-values should be reported for all studies where “tendencies” are claimed, and conclusions should not be drawn from those. Can the authors specify how they decide on reporting non-significant – with no effect VS non-significant – with a tendency? Many results are also based on changes that are close to not meeting significance. Although we understand this is common for vivo work, attention should be paid to the conclusions with high variability, especially based on 3 data points.\n4- In Figure 2, it would be interesting to show the BC/LC ratio as shown in Fig.1 and 3. Although the authors report an increase of Gfp-LC cells in K14-RANK mice compared to K14-Ctl mice (suggesting gain of luminal gene expression markers by BCs and the shift to a luminal identity), such ratio would inform on how the lineages are affected upon enhanced RANK expression in the context of pregnancy.\n5- Can the authors confirm what is used to specify BC vs LC in their analyses (and to specify B/L hybrid)? Is it based on CD24 and CD49 surface expression, K14-K8, etc? Could the change in marker (BC to LC) occur without necessarily affecting a complete switch in lineage?\n6- In figure 3 and supp.3, most panels only show MPA-treated mice and don’t show the untreated condition in parallel. Also, the K14-Ctl contrition is not always shown (e.g., 3c-d-e). These controls would help appreciate the effect of activating RANK on the conclusions drawn (e.g., Ck5, Ck6a, Sma, etc.). For example, it is difficult to appreciate the loss of Ck5 since no condition where Ck5 is expressed is shown.\n7- In Figure 3i, it is important to define how the authors associate genes with chromatin accessibility. Are they looking at proximal promoters? distal enhancers? What rules were considered to define “gene accessibility”?\n8- Although changes in chromatin accessibility are observed between the different models (KANK vs Ctl) and changes in the relative levels of chromatin regulators (eg. Kdm4, Asxl2), and TF accessibility are shown, these observations only allow to draw correlative suggestions. Mechanistic conclusions about epigenetic reprogramming or transcription factor binding/role in lineage specification and tumor development should not be drawn without functional assays allowing to test this (for eg. they have not shown that loss of ASXL2 leads to epigenetic reprogramming and to a loss of basal cell identity and generation of the pre-invasive lesion; other examples also).\n9- In Figure 4f-h, has it been validated that RANK signaling exists in the Gfp+-Control tumors? “These Gfp+ control tumors included Er+ adenocarcinomas and Er- adeno-squamous carcinomas (Fig. 4f-h), indicating that BCs with physiological levels of Rank can generate luminal mammary tumors”. It is not shown that the generation of tumors in this model is due to the physiological level of RANK signaling.\n10- In Figure 4, it is difficult to conclude that the tumor formation is due to the loss of BC lineage fidelity upon enhanced RANK signaling without functionally demonstrating it. Enhanced rank signaling may also impact various pathways that contribute to pro-tumorigenic effects (other than lineage specification).\nMinor comments:\n1- It is unclear what the authors mean by “Sample identity was confirmed by heatmap analysis of known luminal and basal markers” (Supplementary Fig. 2c). Can this be clarified?\n2- The statistical test for panel 2h was not reported.\n3- In Figure 3g-h, the terms accessible and non-accessible should be preferred to open and close chromatin.\n4- At line 272, the sentence should read: “Indeed, one of the most highly accessible chromatin peaks in K14-Rank BCs compared to control BCs corresponded to the promoter of Ehf gene (Fig. 3j), whose mRNA was also found upregulated in the transcriptome of K14-Rank BCs”. The regions/parameters assigned the closest or corresponding genes should be specified.\n5- Can the authors justify using phosphor proteomic analyses to investigate chromatin organization?\n6- Figure 5 is missing panel c to f.\n7- In all GSEA graphs (e.g. 5f), the NES and p-value should be indicated.", "reviewer_3": "This manuscript investigates the impact of overexpressing RANK in basal mammary epithelial cells through a K14 promoter, demonstrating impacts on mammary gland development lactation, and tumorigenesis. Using extensive transcriptomics, proteomics, phospho-proteomics, and ATAC-seq, the authors show that basal-specific RANK creates basal/luminal hybrid states, link these to premalignant lesions and ultimately invasive adenocarcinomas. From their transcriptomic studies, they identify a “basal RANK” signature that correlates with disease outcomes in human breast cancer. The experiments are clearly performed and the data is clear and well presented.\nMajor strengths:\n• Inducible, basal-specific model and the inclusion of fate-mapping elements to track basal-derived cells\n• Comparison with complementary K14-driven RANK ko model provides additional specificy of the observed phenotypes\n• Robust multi-omics approach\n• Generation of DCIS-like preinvasive lesions that can progress to both basal-like and ER+ luminal like cancers\nMajor weaknesses\n• Needs more explanation of rationale for experiments. RANK has been extensively studied in luminal cells as these are the primary producers of this molecule. What is the physiological relevance of basal overexpression of RANK? Was there reason to think that there would be lineage hybridization? It is difficult to assess whether the outcomes followed expectations from previously published literature.\n• While the proteomics data are used to implicate ASXL2 and KDM4A as tumor suppressive epigenetic regulators, and this is relevant to the exploratory studies in human datasets. Without rescue experiments to reconstitute normal ASXL2 or KDM4A levels, causation cannot be assumed.\n• Translational conclusions could be stengthened. Would measuring RANK in human breast cancers lead to identification of metaplasia/lineage hybridization, or improvement of therapy decisions (or alternatively, correlation with established clinical parameters)?", "reviewer_4": "In the manuscript ‘Rank signaling drives basal cell-lineage infidelity leading to mammary tumorigenesis’, by Redondo-Pedraza et al. the authors have carried out a proteomic experiment to compare K14-Rank and K14-Ctrl eight-week-old mammary glands.\nI have reviewed this experiment and the interpretation of the resulting data. My main concern is that the authors have not provided sufficient experimental, analytical and statistical detail, for me to be able to review their data appropriately. As such I am unable to recommend this part of their work for publication.\nFirst, the objective of the experiment is not clearly outlined. Why has it been done in mammary glands, as compared to the luminal and basal cells? I would recommend the authors to improve the clarity of the system employed and specify how this will be helpful in the context of the study.\nIn the methods, we learn that epithelial organoids were harvested for proteomics, and this is followed by a description of the proteomics sample preparation. However, it is not clear how the organoids were treated. The design of the experiment is missing and this needs to be provided as from the main text we learn that you are comparing Rank versus controls.\nFurther, in the Data Analysis section there are key details missing. The authors have used TMTpro labelling, however there is no mention of searching the data with the TMT tags, nor what residues were searched for as fixed or variable modifications. What was the labelling efficiency, as well as partial and overlabelling rates?\nThe authors claim that they have uploaded their data to the PRIDE repository as they report dataset identifiers PXD054102 (proteomics) and PXD058134 (phosphoproteomics). However, they don’t provide reviewer login details. As such these data cannot be reviewed and this is a major issue as it does not comply with the requirements of Data availability of Nature Communications.\nIn general this experiment demonstrated weak quantitative trends on the protein and phosphoprotein level, both in terms of effect size and significance. The authors put forward a fold change cutoff of +/-1.2 which is very low indeed. There wouldn’t be anything problematic with small effect sizes if the authors would be able to demonstrate that they have sufficient statistical power to detect these small fold changes. The authors ought to provide evidence on this aspect given the technical variability of their system and analytical measurement. Moreover, details of replication are missing from the Methods and needs to be specified. In terms of significance, are the authors reporting raw or adjusted p-values – again I am unable to review their FDR considerations. Looking at Supplementary table 4, I could not find the quantification of the authors’ phosphoproteomics data. In the protein quantification table, there is Reporter_intensity_1_TotalProt (10 columns) and then again Reporter_intensity_1_TotalProt (18 columns). This is misleading as in the Methods the authors write about using TMT pro 18 plex. In the final statistical tables, I can see the experimental Contaminats such as keratins within their statistically significant list of proteins and phosphosites. I would recommend the authors to filter these out from the final lists of candidates, as they are clearly not candidates but experimental contaminants. Further details that are required to be provided before this experiment can be considered, are phosphosite localisation probabilities and apparent motif annotation that has been employed.\nMinor comments\nMouse protein database (UniProtKB/TrEMBL; 21,990 sequences), needs to be dated."}, "author_response": "We thank the reviewer for this question. Results from our group and others have demonstrated that RANKL expression in the MG is regulated by progesterone and it is restricted to hormone sensing positive cells, which express hormone receptors (progesterone and estrogen receptors) (reviewed in (Slepicka et al., 2021; Perez-Chacon et al., 2025)) and indicated in the manuscript (lines 71-73). Indeed, analyses of published scRNA-seq datasets confirmed that the natural source of RANKL within the epithelial compartment of the human breast and mouse mammary gland is the hormone-sensing luminal cell (LC) population and this is also true for BRCA1 mutation carriers (TNFSF11 Fig. R1.1a) (Bach et al., 2017; Reed et al., 2024). This has now been indicated in the revised manuscript (Supplementary Fig. 2i, j in the revised manuscript and lines 158-159). Additionally, IHC staining of K14-Ctrl/K14-Rank MGs confirmed that Rankl protein expression is restricted to the LC compartment upon MPA treatment (Supplementary Fig. 3f in the revised manuscript, lines 206-207).\nIn human ductal carcinomas in situ (DCIS) and some early-stage luminal A breast adenocarcinomas, RANKL expression is still restricted to HR+ luminal cells. We confirmed in published scRNA-seq datasets that the primary source of RANKL in DCIS is the luminal compartment and the immune infiltrating cells and that its expression was reduced in invasive ductal carcinomas (IDC) (Fig. R1.1b) (Tokura et al., 2022). Moreover, we analysed the expression of RANKL in a separated cohort of 57 DCIS and 313 IDC samples (Bergholtz et al., 2020) (Supplementary Fig. 10a, b in the revised manuscript, lines 385-391).\nFigure R1.1. (a) RANKL (TNFSF11) expression profile obtained from human BRCA1 carriers scRNA-seq databases (Reed et al., 2024). BC: basal cells; LASP: luminal associated secretory progenitors; LHS: luminal hormone sensing. (b) RANKL (TNFSF11) expression profile obtained from human DICS and IDC scRNA-seq databases (Tokura et al., 2022). Luminal cells and immune cells are the primary source of RANKL in human DICS.\n2. The drastic changes in chromatin accessibility upon enhanced Rank signaling is a very interesting finding. While the authors demonstrated these chromatin remodeling changes are associated with various transcription factors and epigenetic regulators like Asxl2, they should verify which of these are the critical regulator(s) of this chromatin reorganization. This can be validated in vitro using mouse/human mammary or breast cancer cell lines through targeting these transcription factors and epigenetic regulators under enhanced Rank signaling.\nWe agree with the reviewer that underpinning which are the crucial epigenetic regulators that operate downstream of RANK and are ultimately responsible for the chromatin reorganization would be highly relevant. We have tried this in different cellular models where we could manipulate the activity of RANK pathway and the expression of epigenetic regulators.\nWe have mainly focused our efforts on mouse BC44 and human MCF10a cells as plausible models of basal non-transformed basal cells.\nHowever, in these cell lines Rank overexpression did not induce consistent changes in B/L identity. We also determined that in these cellular systems Rank does not directly reduce the expression of Asxl2, Kmd4 nor increases Batf3 (Fig. R1.3a, b).\nFigure R1.2. (a) Gene expression analysis of Rank (Tnfrsf11a), Krt14, Ehf, Krt8, and Krt5 in BC44 cells with ectopic Rank expression compared to empty vector (EV). (b) Gene expression analysis of RANK (TNFRSF11A), KRT14, EHF, KRT8 and KRT5 in MCF10A with ectopic Rank expression compared to empty vector (EV) controls. Each dot represents a technical replicate. One representative experiment of 3 (a) and 4 (b) is shown. Mean, SEM and p-values are shown (*, p<0.05; **, 0.001<p<0.01; ***, 0.0001<p<0.001; ****, p<0.0001).\nAlso, as the referee suggested, we overexpressed Asxl2 and Kdm4a in BC44 cells and assessed the impact of Rank overexpression in this context. We achieved overexpression of both genes as well as Rank (Fig. R1.4a). No consistent gene expression changes in basal and luminal markers were found (Fig. R1.4a). However, immunofluorescence analyses revealed a trend of increased Ck5 expression in BC44 Rank overexpressing cells with co-expression of Asxl2 and Kdm4a suggestive of increased basal cell identity in line with in vivo mouse observations (Fig. R1.4b).\nWe next assessed the impact of the downregulation of Asxl2 and Kdm4a by shRNA. Downregulation of Kdm4a severely impacted BC44 cell viability and proliferation, therefore we could not pursue additional studies in this context. Asxl2 silencing was modestly achieved (Fig. R1.5a) and in this case, Krt5 and Krt14 were downregulated in Rank-overexpressing conditions, while Ehf and Krt8 expression increased (Fig. R1.5a). However, minimal differences in Ck5 protein levels were found between conditions as determined by immunofluorescence analyses (Fig. R1.5b).\nSince none of the previous experiments provided clear results regarding changes in cell lineage identity, we then tried to assess whether targeting the epigenetic machinery would alter Rank-induced B/L transition. We reasoned that regulation of H3K27me3 by Rank signaling may alter lineage specification and B/L conversion given that global levels of H3K27me3 accompany mammary lineage commitment (Holliday et al., 2018). Therefore, we treated BC44 cells with the epigenetic drugs GSK-126 and GSK-J4 and analyzed changes in cell identity genes. GSK-126 suppresses H3K27me3 by inhibiting EZH2 methyltransferase, a member of the PRC2 complex, that promotes H3K27me3, while GSK-J4 hinders H3K27 demethylation by selectively inhibiting the KDM6 histone demethylase family (Heinemann et al., 2014). We assessed changes in lineage specification by analyzing the expression of basal (Krt5/Krt14) and luminal markers (Krt8) or LC-associated gene Ehf by qPCR. In none of the treatments, were we able to observe consistent and reproducible changes in the lineage-associated markers (Fig. R1.6).\nAltogether these results demonstrate that these in vitro cell models are not appropriate to assess the mechanism underlying the B/L conversion driven by Rank in vivo, and therefore we have not included them in the revised manuscript.\nFigure R1.6. (a-b) Gene expression analysis of Rank (Tnfrsf11a) Asxl2, Kdm4a, Krt5, Krt14, Krt8 and Ehf BC44 cells with ectopic Rank expression compared to empty vector (EV) treated with GSK-126 (a) and GSK-J4 (b) Each dot represents a technical replicate. One representative experiment out of 3 is shown. Mean, SEM and significant p-values for Rank overexpressing cells treated or not with inhibitors are shown (*, p<0.05; **, 0.001<p<0.01; ***, 0.0001<p<0.001; ****, p<0.0001).\nConsidering that global levels of the methylation mark H3K27me3 accompany mammary lineage commitment (Holliday et al., 2018), we decided to address the status of the H3K27me3 mark in our models and unveiled that hybrid K14-Rank cells show levels of H3K27me3 intermediate to those of luminal and basal cells (Supplementary Fig. 9a the revised manuscript). We then reasoned that Rank-driven regulation of H3K27me3 may contribute to B/L conversion and the appearance of the hybrid B/L cells in vivo. To test this hypothesis, we treated the animals with H3K27me3 modulators that promote a rewiring of the chromatin accessibility (Supplementary Fig. 9b in the revised manuscript).\nBriefly, morphological assessment and flow cytometry analyses revealed that GSK-J4 treatment increased the frequency of Gfp+ LCs (Fig. 4h and Supplementary Fig. 9e in the revised manuscript), suggesting that modulation of H3K27me3 promotes a complete basal to luminal conversion. Importantly, the frequency of Gfp+ B/L hybrid cells was reduced in K14-Rank MGs treated with GSK-J4, while the opposite trend was observed upon GSK-126 treatment (Fig. 4i). Moreover, a drastic reduction of preinvasive lesions containing Ck14+Ck8+ hybrid cells was found in the GSK-J4-treated K14-Rank MGs and most lesions contained lineage-restricted BC Ck14+ Ck8- cells. In K14-Rank MGs treated with GSK-126, a small reduction of lesions containing hybrid cells was also found (Fig. 4j, k in the revised manuscript). These new results confirm that basal Rank signaling favors the onset of B/L hybrid cells by remodeling the epigenome, since targeting chromatin modifiers with drugs impacts this B/L conversion These results are now included as new panels included (Fig. 4g-k and Supplementary Fig. 9) in the results section and discussed in the revised manuscript. Please see results section lines 344-380 and lines 474-483 in the discussion section for a detailed description.\n3. Additionally, the clinical relevance of their findings can be enhanced with further confirmation of whether these chromatin remodeling changes and Rank transcriptional signature can be similarly induced in human invasive or luminal breast cancer cells upon Rankl stimulation.\nResults included in Fig. 5 and Supplementary Fig. 10 from the revised manuscript highlight the clinical relevance of our findings. As preinvasive lesions formed by hybrid cells in K14-Rank mice resemble ductal carcinoma in situ (DCIS), we have assessed expression patterns of Rank-BC transcriptional signature, as well as TF and epigenetic modifiers identified in our mouse models in a valuable cohort of 57 DCIS and 313 IDC samples classified by PAM50 (Bergholtz et al., 2020).\nImportantly Rank-BC signature was higher in IDC compared to DCIS (Fig. 5a, revised manuscript) and, attending to PAM50 classification, show higher levels in HER2 DCIS (Fig. 5b, revised manuscript), which associate with increase recurrence rates (Meyer et al., 2025). Indeed, we have found that basal Rank gene signature is predictive of DCIS recurrence and progression in an independent cohort of DCIS (Fig. 5e in the revised manuscript). As stated in the manuscript this basal Rank signature can be exploited to assess which DCIS is likely to recur, favoring tailored therapies while avoiding unnecessary treatments.\nAttending to the request of the referee, further data mining in the transcriptomic data from the cohort of DCIS and IDC samples (Bergholtz et al., 2020) unveiled an opposite association of EZH2 and KDM6A expression with “basalness”, suggesting that they can modulate B/L identity even in IDC, in line with changes in mammary cell identity found upon treatment with their pharmacological inhibitors, GSK-126 and GSK-J4, in the K14-Rank preclinical model. Moreover, similar to ASXL2, a strong positive association between RANK and EZH2 was found in IDC, opposite to DCIS. These results support our conclusions involving basal RANK in the regulation of B/L cell identity through epigenetic mechanisms in human breast cancer. These results are now included in the main manuscript (lines 396-407) and in Fig. 5c and in Supplementary Fig. 10c.\nNext, per request of the referee, we also analyzed the expression of the basal Rank gene signature (Rank-BC) in transcriptomic data from breast cancer patient-derived orthoxenografts (PDXs) treated in vivo with RANKL or the inhibitor Rank-Fc from our RANK+ PDX cohort (Ciscar et al., 2023). As shown in Fig. R1.7a, despite no significant differences were detected, we observed an opposite trend between RANKL and Rank-Fc treatments, suggesting that the modulation of RANK signaling may slightly impact cell identity even in full blown breast cancer. This is striking considering the multiple chromosomal abnormalities already present in breast cancer cells.\nMoreover, exacerbation of Rank signaling led to significant changes in several pathways associated with chromatin remodeling and epigenetic regulation in both PDX models (untreated vs RANKL-treated) and upon RANK pathway inhibition (untreated vs Rank-Fc) (Fig. R1.7b-d). Although the number of PDXs analyzed is limited, these results suggest that Rank signaling induces chromatin remodeling changes in breast adenocarcinomas, which may modulate tumor cell identity. We decided not to include these results in the revised manuscript, to maintain the focus in breast cancer initiation, as we have demonstrated that chromatin remodeling changes and loss of basal cell identity driven by Rank play a central role in these early stages of transformation. Future studies will be required to properly assess the potential of Rank pathway to induce chromatin remodeling changes in breast adenocarcinomas and its functional and therapeutic consequences.\nCompared to the corresponding controls, the luminal breast cancer cells MCF7 with ectopic RANK expression tended to show lower levels of ASXL2, lower KDM4A and a drastic reduction of BATF3, accompanied by trends pointing towards a reduction in the basal markers KRT5/KRT14 and an increase in the luminal KRT8 and EHF genes, emulating most of the changes observed in vivo. However, RANKL did not have a consistent effect in the regulation of basal/luminal markers.\nSKBr3 cells with RANK ectopic expression show higher levels of ASXL2, KDM4A, BATF3, KRT8 and EHF, than corresponding controls. RANKL stimulation led to a clear increase in the expression of the basal KRT14 and KRT5 and moderate decrease in EHF.\nFinally, in Hs578T cell line the main changes observed were an increase in the expression levels of KDM4A, KRT14, KRT8 and EHF upon RANK overexpression, while RANKL treatment reduced KRT14 and increased EHF levels.\nThe fact that no common patterns were found among breast cancer cell lines upon RANK expression or RANKL treatment and that changes observed were modest, evidenced that these in vitro models cannot recapitulate the in vivo breast cancer complexity and have not been included in the revised manuscript.\n4. Finally, the authors should perform functional studies of the target genes of the key regulator(s) found to underscore these chromatin remodeling changes to explain their biological roles in mediating the basal-luminal cell conversion and tumorigenesis.\nPlease see the response to question 2 above where we detail a range of experiments we performed, and the results obtained to address this point both in vitro and in vivo settings.\nBach K, Pensa S, Grzelak M, Hadfield J, Adams DJ, Marioni JC, Khaled WT. Differentiation dynamics of mammary epithelial cells revealed by single-cell RNA sequencing. Nat Commun 8 2128 (2017). https://doi.org/10.1038/s41467-017-02001-5\nBergholtz H, Lien TG, Swanson DM, Frigessi A; Oslo Breast Cancer Research Consortium (OSBREAC); Daidone MG, Tost J, Wärnberg F, Sørlie T. Contrasting DCIS and invasive breast cancer by subtype suggests basal-like DCIS as distinct lesions. NPJ Breast Cancer. 2020 Jun 17;6:26. doi: 10.1038/s41523-020-0167-x. PMID: 32577501; PMCID: PMC7299965.\nCiscar M, Trinidad EM, Perez-Chacon G, Alsaleem M, Jimenez M, Jimenez-Santos MJ, Perez-Montoyo H, Sanz-Moreno A, Vethencourt A, Toss M, Petit A, Soler-Monso MT, Lopez V, Gomez-Miragaya J, Gomez-Aleza C, Dobrolecki LE, Lewis MT, Bruna A, Mouron S, Quintela-Fandino M, Al-Shahrour F, Martinez-Aranda A, Sierra A, Green AR, Rakha E, Gonzalez-Suarez E. RANK is a poor prognosis marker and a therapeutic target in ER-negative postmenopausal breast cancer. EMBO Mol Med. 2023 Apr 11;15(4):e16715. doi: 10.15252/emmm.202216715. Epub 2023 Mar 7. PMID: 36880458; PMCID: PMC10086586.\nHeinemann B, Nielsen JM, Hudlebusch HR, Lees MJ, Larsen DV, Boesen T, Labelle M, Gerlach LO, Birk P, Helin K. Inhibition of demethylases by GSK-J1/J4. Nature. 2014 Oct 2;514(7520):E1-2. doi: 10.1038/nature13688. PMID: 25279926.\nHolliday H, Baker LA, Junankar SR, Clark SJ, Swarbrick A. Epigenomics of mammary gland development. Breast Cancer Res. 2018 Sep 3;20(1):100. doi: 10.1186/s13058-018-1031-x. PMID: 30176939; PMCID: PMC6122685.\nMeyer L, Jackson HW, Eling N, Zhao S, Usui G, Dakhli H, Schraml P, Dettwiler S, Elfgen C, Varga Z, Moch H, de Souza N, Bodenmiller B. A stratification system for breast cancer based on basoluminal tumor cells and spatial tumor architecture. Cancer Cell. 2025 Sep 8;43(9):1637-1655.e9. doi: 10.1016/j.ccell.2025.06.019. Epub 2025 Jul 17. PMID: 40680742.\nPalafox M, Ferrer I, Pellegrini P, Vila S, Hernandez-Ortega S, Urruticoechea A, Climent F, Soler MT, Muñoz P, Viñals F, Tometsko M, Branstetter D, Dougall WC, González-Suárez E. RANK induces epithelial-mesenchymal transition and stemness in human mammary epithelial cells and promotes tumorigenesis and metastasis. Cancer Res. 2012 Jun 1;72(11):2879-88. doi: 10.1158/0008-5472.CAN-12-0044. Epub 2012 Apr 10. PMID: 22496457.\nPérez-Chacón G, Santamaría PG, Redondo-Pedraza J, González-Suárez E. RANK/RANKL Signaling Pathway in Breast Development and Cancer. Adv Exp Med Biol. 2025;1464:309-345. doi: 10.1007/978-3-031-70875-6_16. PMID: 39821032.\nReed, A.D., Pensa, S., Steif, A. et al. A single-cell atlas enables mapping of homeostatic cellular shifts in the adult human breast. Nat Genet 56, 652–662 (2024). https://doi.org/10.1038/s41588-024-01688-9\nSlepicka, P. F., Somasundara, A. V. H., & Dos Santos, C. O. (2021). The molecular basis of mammary gland development and epithelial differentiation. Semin Cell Dev Biol, 114, 93-112. https://doi.org/10.1016/j.semcdb.2020.09.014\nTokura M, Nakayama J, Prieto-Vila M, Shiino S, Yoshida M, Yamamoto T, Watanabe N, Takayama S, Suzuki Y, Okamoto K, Ochiya T, Kohno T, Yatabe Y, Suto A, Yamamoto Y. Single-Cell Transcriptome Profiling Reveals Intratumoral Heterogeneity and Molecular Features of Ductal Carcinoma In Situ. Cancer Res. 2022 Sep 16;82(18):3236-3248. doi: 10.1158/0008-5472.CAN-22-0090. PMID: 35852797.\n14"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have adequately addressed my comments in the revised manuscript.", "reviewer_2": "The authors have adequately addressed my comments. Modifications and new data have been incorporated into the revised version, addressing particular concerns regarding mechanistic and functional validations.", "reviewer_3": "Thank you for the thorough and transparent revision. In my view the authors have substantially addressed the three major concerns from the first round.\n1. Physiological rationale for basal RANK and expectation of lineage hybridization: The revised text clarifies prior evidence for higher RANK in basal cells and for RANK dependent basal bipotency. Importantly, the new basal RANK loss of function carcinogenesis experiment (K14 RankΔ/Δ + MPA/DMBA) shows lower tumor incidence, delayed onset, and GFP negative tumors, arguing that physiological basal RANK is required for basal derived tumorigenesis. This directly addresses the request for physiological relevance.\n2. Causality of ASXL2/KDM4A vs. chromatin state mechanism: The authors report extensive in vitro rescue attempts (BC44/MCF10A; ASXL2/KDM4A overexpression/knockdown with RANK OE) that were modest or inconsistent, concluding these systems do not recapitulate the in vivo B/L switch. They then use in vivo epigenetic modulation to show hybrid cells have intermediate H3K27me3 and that GSK J4 (↑H3K27me3) increases GFP⁺ luminal cells and reduces GFP⁺ hybrid cells and hybrid cell lesions, with GSK 126 showing the opposite trend. These data mechanistically strengthen the chromatin state model even though a gene specific in vivo rescue for ASXL2/KDM4A remains outstanding. I recommend tempering any language that implies those two factors are the causal nodes and instead framing them as leading candidates within an H3K27-centered mechanism.\n3. Translational significance: The revision adds human and PDX analyses: the Rank BC (basal RANK) signature is higher in IDC vs. DCIS, enriched in HER2 DCIS, and predicts DCIS recurrence in an independent cohort; in IDC, RANK–EZH2 association and opposing EZH2/KDM6A–“basalness” patterns support an epigenetic axis; PDXs show directionally consistent effects with RANKL vs. RANK Fc. These additions strengthen the translational narrative, with appropriate acknowledgement that prospective validation will be needed before clinical use.\nMinor suggestions:\n1. Where discussing ASXL2/KDM4A, adjust phrasing from “causal regulators” to “candidate mediators” and point explicitly to the in vivo H3K27me3 manipulation as the strongest mechanistic evidence currently available.\n2. Keep a one sentence limitation noting that available in vitro systems did not reproduce the lineage switch, explaining the pivot to in vivo pharmacologic modulation.", "reviewer_4": "The authors have addressed my concerns with regards to methodological description. A few remaining points that should be addressed before I can recommend acceptance:\n• The rebuttal claims that the proteomic experiments have been detailed in line 331-337. Unfortunately, I cannot identify changes in the text that would point to improved clarity of the proteomic experiment.\n• Unfortunately, I have identified flaws in the authors’ response with regards to their interpretation of the quantification and statistical analysis used. It is incorrect to adjust fold-changes based on the weaknesses of the TMT-based MS2-level quantification strategy. The highlighted study by Martinez-Val et al., 2016 clearly describes the problem. In order for the authors to be able to select arbitrary cutoffs, especially ones that allow very small effect sizes, they would need to provide evidence that they have the statistical power to detect such minimal effect sizes. At no point in the Martinez-Val et al study, was it recommended to lower fold-change thresholds or to use specific cutoffs. It is wrong to implement cutoffs from studies employing different experimental systems with different technical and biological variabilities. On the contrary, they recommend: ‘On the other hand, for ratios smaller than 2.5, a prior analysis of the statistical power could be done. This will be useful to determine the appropriate number of replicates.’ Hence, the author's interpretation of the Martinez-Val et al study and the arguments used to choose the arbitrary low cutoffs ±1.2 fold change (log2 Fc of -/+ 0.3) are not sound . I recommend the authors to employ an analysis that would identify trends in pathways rather than relying on individual proteins with such small fold changes.\n• It is unclear where the following has been added: ‘These findings confirm that our TMT labelling efficiency was high, with minimal unintended modifications at STY sites, ensuring the robustness of our quantitative proteomics workflow. We will provide a more detailed breakdown of these percentages and their implications in the revised manuscript.’\n• Please remove the following duplication:\n‘750 independent biological replicate, derived from distinct mice. After labeling, reactions\n751 were quenched with 5% hydroxylamine. The labeling reaction was quenched by\n752 adding 5% hydroxylamine.’"}, "author_response": "We have followed the reviewer suggestion and modified accordingly some sentences in the Discussion section (lines 481-490). Changes are now indicated in green.\n2. Keep a one sentence limitation noting that available in vitro systems did not reproduce the lineage switch, explaining the pivot to in vivo pharmacologic modulation.\nWe have tried to incorporate this suggestion in our manuscript. However, including the sentence “available in vitro systems did not reproduce the lineage switch” would have required to indicate “data not shown” next to it, which we do not find suitable due to current publication policies.\nTherefore, we have reinforced in the discussion the in vivo pharmacologic modulation (lines 482-490) that we hope the reviewer will find appropriate.\nWe want to thank the reviewer for the comprehensive and fruitful review performed. We hope that we have now addressed the minor concerns raised in a convincing manner.\nWe believe that our manuscript has improved greatly and our results are now more relevant thanks to the suggestions made by the reviewers in general and reviewer n.3 in particular.\nReviewer #4 (Remarks to the Author):\nThe authors have addressed my concerns with regards to methodological description. A few remaining points that should be addressed before I can recommend acceptance:\nWe sincerely thank Reviewer #4 for the careful evaluation of our manuscript and for the constructive comments that significantly improved the clarity and rigor of the proteomic analyses. We have now fully addressed all remaining concerns, as detailed below.\n• The rebuttal claims that the proteomic experiments have been detailed in line 331-337. Unfortunately, I cannot identify changes in the text that would point to improved clarity of the proteomic experiment.\nWe thank the reviewer for pointing out that the previous version did not clearly reflect the proteomic workflow in the main text. In the revised manuscript, we have now added explicit and detailed clarification in the Results and Methods section. In particular, lines 331-346 in the Results section now read: “To investigate how enhanced basal Rank signaling remodels epithelial protein programs, we performed quantitative TMT-based proteomic and phosphoproteomic analyses on freshly isolated epithelial organoids from virgin 8-week-old K14-Rank and K14-Ctrl females, which reflect the Rank gain-of-function epithelium”\nIn lines 746-755 and 767-781 of the methods section additional clarifications have been included explaining the rationale, sample source, and experimental design of the proteomic workflow.\nThis addition should ensure that readers immediately understand the biological rationale and sample design underlying the proteomic experiment, directly addressing the reviewer’s concern about clarity.\nFurthermore, the Methods section has been substantially expanded to provide full technical detail of the proteomic workflow. The full description now appears in:\n“Sample preparation for proteomic analysis” (lines 747-765),\n“TMTpro Labeling” (lines 767-781),\n“Data Analysis” (lines 819-855).\nTogether, the main text and Methods now provide a transparent, complete, and traceable description of the proteomic and phosphoproteomic experiments.\n• Unfortunately, I have identified flaws in the authors’ response with regards to their interpretation of the quantification and statistical analysis used. It is incorrect to adjust fold-changes based on the weaknesses of the TMT-based MS2-level quantification strategy. The highlighted study by Martinez-Val et al., 2016 clearly describes the problem. In order for the authors to be able to select arbitrary cutoffs, especially ones that allow very small effect sizes, they would need to provide evidence that they have the statistical power to detect such minimal effect sizes. At no point in the Martinez-Val et al study, was it recommended to lower fold-change thresholds or to use specific cutoffs. It is wrong to implement cutoffs from studies employing different experimental systems with different technical and biological variabilities. On the contrary, they recommend: ‘On the other hand, for ratios smaller than 2.5, a prior analysis of the statistical power could be done. This will be useful to determine the appropriate number of replicates.’ Hence, the author's interpretation of the Martinez-Val et al study and the arguments used to choose the arbitrary low cutoffs ±1.2 fold change (log2 Fc of -/+ 0.3) are not sound . I recommend the authors to employ an analysis that would identify trends in pathways rather than relying on individual proteins with such small fold changes.\nWe agree with the reviewer’s concerns about the interpretation of MS2-level TMT quantification. In the revised manuscript, we have completely eliminated the use of fold-change cutoffs (including the ±1.2 or log₂ ±0.3 threshold) and no longer use p-values or multiple-testing–adjusted thresholds to define differentially expressed proteins.\nImportantly, we have also updated Supplementary Table 4 to remove any columns or annotations implying differential expression, statistical significance, or FC-based filtering. It now contains only the unfiltered quantitative matrices and ranked values, consistent with the revised analytical strategy.\nIn response to the reviewer’s recommendation, we now rely exclusively on rank-based pathway analysis to capture coordinated biological trends rather than attempting to call individual differentially expressed proteins. For the proteome, proteins are ranked by their log₂ fold change (K14-Rank vs K14-Ctrl) and analyzed using GSEA with the full quantified proteome as background.\nFor the phosphoproteome, we have replaced the previous site-level significance filtering with a gene-level ranking approach: quantified phosphosites are collapsed to a single value per gene (maximum absolute log₂ fold change), genes are ranked by this signed fold change, and the resulting list is analyzed by GSEA using GO Biological Process and Reactome collections. This preserves the reviewer’s requirement to avoid arbitrary thresholds while still allowing biologically interpretable pathway-level analysis.\nThe analysis has been described in Methods (line 819-855).\nTogether, these changes constitute a complete restructuring of the quantitative analysis framework, fully aligned with the reviewer’s concerns and the recommendations of Martinez-Val et al. (2016), and avoid any overinterpretation of small effect sizes inherent to MS2-level TMT data.\n• It is unclear where the following has been added: ‘These findings confirm that our TMT labelling efficiency was high, with minimal unintended modifications at STY sites, ensuring the robustness of our quantitative proteomics workflow. We will provide a more detailed breakdown of these percentages and their implications in the revised manuscript.’\nWe thank the reviewer for pointing out that, in the previous revision, the information regarding TMTpro labelling efficiency and STY over-labelling was not clearly identifiable in the manuscript. We have now explicitly incorporated this information into the Data Analysis section of the Methods (lines 819-855) in the revised manuscript). In this new text, we report that TMTpro labelling efficiency exceeded 98%, that <1% of total precursor intensity corresponded to unlabelled peptides, and that over-labelling at serine, threonine, and tyrosine residues accounted for 4.99% of precursor intensity (6.46% of PSMs), values that are lower than typically reported for TMTpro datasets. We also added the conclusion that these metrics confirm high labelling efficiency with minimal unintended STY modifications, ensuring the robustness of our quantitative proteomics workflow. These details now appear clearly in the revised Methods and directly address the reviewer’s concern.\n• Please remove the following duplication:\n‘750 independent biological replicate, derived from distinct mice. After labeling, reactions\n751 were quenched with 5% hydroxylamine. The labeling reaction was quenched by\n752 adding 5% hydroxylamine.’\nWe thank the reviewer for noting the duplicated description of the hydroxylamine quenching step. This duplication has now been fully removed. The Methods section (“TMTpro Labeling”, line 771) contains a single, correct sentence describing quenching with 5% hydroxylamine, and no redundant text remains."}, {"round": 3, "reviewer_comments": {"reviewer_3": "The authors have addressed my concerns fully.", "reviewer_4": "The authors have improved the clarity of their experiment and modified the analysis of their proteomic and phosphoproteomic data by employing ‘no fixed fold-change or significance thresholds to define differentially expressed proteins’. Keeping in line with this, they ought to remove their previous statistical analysis from Fig 4F as it is no longer described in the methods and goes against the figure legend description (i.e. no cutoffs should be highlighted in red/blue and no p-values should be presented). I would suggest to use a simple waterfall plot of their ranked fold changes where they can highlight protein candidates.\nIn terms of their GSEA analysis presented in Supplementary Table 4, it is unclear whether this relates to Fig 4F or Supplementary Fig 8G? The methodology about how the proteomic and phosphoproteomic data was used to create the GSEA output needs to be added and should include details of how the ranking was done for phosphoproteins that were also detected at the protein level. There is a mention in the text about applying ‘GSEA to the full ranked protein list and phosphosite rankings’, however, the legend of Supplementary Fig 8G only refers to ‘GSEA of phosphoproteome changes’."}, "author_response": "As suggested by the reviewer, the previous statistical analysis has been removed from Fig. 4F.\nFig. 4F no longer shows highlighted candidates based on cutoffs and p-values.\nIn terms of their GSEA analysis presented in Supplementary Table 4, it is unclear whether this relates to Fig 4F or Supplementary Fig 8G? The methodology about how the proteomic and phosphoproteomic data was used to create the GSEA output needs to be added and should include details of how the ranking was done for phosphoproteins that were also detected at the protein level. There is a mention in the text about applying ‘GSEA to the full ranked protein list and phosphosite rankings’, however, the legend of Supplementary Fig 8G only refers to ‘GSEA of phosphoproteome changes’.\nTable 4 compiles results for all the proteomic/phosphoproteomic analyses, as indicated in the “readme” and the legends of every sheet. GSEA has been applied to to the full ranked protein list and phosphosite rankings’. Supplementary Fig. 8G represents GSEA of phosphoproteome changes’ as indicated.\nFor proteome-level GSEA, all quantified proteins were ranked solely by their log₂ fold change (K14-Rank vs K14-Ctrl), without applying fold-change cutoffs or statistical significance thresholds, and the full ranked list was used as input for GSEA.\nFor phosphoproteome-level GSEA, the analysis was performed independently of protein-level abundance changes. Specifically, phosphosite-level quantifications were first collapsed to a single value per gene to enable gene-centric pathway analysis. For each gene, the phosphosite exhibiting the maximum absolute log₂ fold change was selected, and the signed log₂ fold change of that site was used to rank genes. This approach captures the dominant direction and magnitude of phosphorylation change per protein while avoiding averaging across sites that may change in opposite directions.\nImportantly, protein-level fold changes were not used to weight, normalize, or adjust phosphoproteomic rankings, even when the corresponding protein was also detected in the global proteome dataset. Proteome and phosphoproteome GSEA analyses were therefore conducted as conceptually independent analyses, reflecting changes in protein abundance and phosphorylation state, respectively.\nAll GSEA analyses were performed using the full ranked lists and appropriate background sets, consistent with a rank-based framework designed to identify coordinated pathway-level trends rather than individual differentially expressed entities. These methodological details have now been added explicitly to the Methods section Proteomics and Phosphoproteomics Data Analysis to ensure full transparency and reproducibility (line 820-868)."}]

Decision letter