A hormetic transcriptional program coregulates invasion, proliferation and dormancy to define metastatic potential

Domain: chemistry
Published: 2026-03-04
Editorial outcome: major_revision

Abstract
Metastasis remains the leading cause of cancer mortality, yet the factors that determine metastatic competence and the timing of its acquisition are not well understood. Using human breast cancer cohorts and multi-omics in preclinical models, here we show that cells exit the primary tumor with pre-set metastatic potential. Integrating spatial transcriptomics with single-cell gene expression and chromatin profiles, we identify the transcription factor Prrx1 as a master regulator of dissemination: beyond its known role in invasion, it represses proliferation by acting on cyclins and cell-cycle inhibitors (Ccnd1/2, Cdkn2a/b/c) and activates a dormancy program (Gas6, Mme, Ogn). Intermediate Prrx1 levels optimize the trade-off between invasion and proliferation/dormancy, producing a hormetic (nonlinear) relationship between Prrx1 expression and metastatic burden. Combined invasion and proliferation signatures strongly stratify breast cancer prognosis, underscoring the clinical relevance of this phenotypic interplay. These findings indicate that metastatic competence emerges from the integration of Prrx1 co-regulated molecular programs, explaining why some invasive cells enter dormancy while others resume growth.

Introduction
Despite significant advances in cancer treatment, which have extended median survival rates, most patients ultimately succumb to metastatic relapse1. Metastasis is a highly complex, multistep process that only a small fraction of cancer cells manages to complete2. The timing of metastasis is highly variable, some patients are diagnosed with metastatic disease at the time of detection, while in others, metastases appear years after initial treatment3–5 . Recent advances have shed light on the biology of metastatic cells, highlighting cell-autonomous characteristics and microenvironmental effectors as key determinants of metastatic success6. Additionally, there is little evidence supporting a strong genetic basis for metastasis7–9 suggesting that cancer cell heterogeneity and plasticity may underlie metastatic behavior10. However, the biological features that define metastatic cells and when they are acquired remain poorly understood.

One possibility is that pre-existing cell states influence the future behavior of metastatic cells, due to pure intrinsic characteristics or distinct responsiveness to environmental signals at the metastatic site11. Epithelial-mesenchymal transition (EMT) is a well-known driver of cell plasticity and tumor dissemination12,13. EMT is a dynamic process in which cells can lie at a spectrum of intermediate states between the epithelial and mesenchymal phenotypes14. Different EMT states in cancer cells have been associated with varying metastatic potential, with hybrid E/M cells often exhibiting the highest metastatic competence15–18.

Recently, we described an EMT-driven intratumor heterogeneity in breast cancer, identifying two distinct trajectories, respectively associated with invasion or inflammation, each controlled by different EMT transcription factors19. We found that among them, Prrx1 is required for invasion and consequently, metastasis. However, Prrx1 high expression prevents metastatic colonization20,21. These contrasting results suggest that the integration of different cellular programs that shape metastatic potential may be influenced by Prrx1 levels.

By integrating patient data, single cell sequencing modalities, genetic mouse models bearing tumors with distinct levels of Prrx1 and engineered cell lines, we demonstrate that Prrx1 expression in the primary tumor dictates metastatic competence by simultaneously regulating key cellular traits: invasion, proliferation, and dormancy. In addition, we found that Prrx1 levels exhibit a hormetic behavior, a biphasic dose–response in which intermediate levels promote metastasis, while high levels have an inhibitory effect22,23, reconciling the apparently conflicting data mentioned above19–21. Furthermore, we identify Prrx1 as a master regulator that integrates invasion, proliferation and dormancy programs and show that a combinatorial approach using molecular signatures derived from our analyses enables the identification of a breast cancer subgroup with poor prognosis.

Results
Non-linear correlation between PRRX1 levels and metastatic incidence in breast cancer patients

We had previously shown that high levels of Prrx1 expression in the primary tumor correlated with the lack of metastasis in breast cancer patients20. To better assess the relationship between PRRX1 expression in cancer cells and metastatic incidence without the influence of stromal derived PRRX1 expression, we analyzed 94 (grade 3) infiltrating ductal carcinoma (IDC) tumors in an adjuvant setting from patients who underwent surgery between 2010 and 2011. PRRX1 was visualized using immunohistochemistry (INHC) in combination with an anti-keratin antibody cocktail to specifically identify cancer cells, and intratumor heterogeneity and intensity of PRRX1 signal in cancer cells was evaluated by a pathologist. Based on both the intensity and heterogeneity of PRRX1 expression in cancer cells, tumors were classified into three groups: the majority (82%, with no detectable PRRX1-positive cancer cells) as PRRX1-negative; 10% as PRRX1-intermediate, characterized by a sparse presence of cells expressing intermediate PRRX1 levels; and the remaining 8% as PRRX1-high, containing a substantial number of cancer cells strongly positive for PRRX1 (Fig. 1a). To substantiate the existence of these three groups in a quantitative manner, a subset of samples containing all groups were analyzed using image-analysis software on a single-cell basis (see Methods). Unsupervised clustering of the patients’ scores based on (i) percentage of PRRX1+ and (ii) intensity of PRRX1 in cancer cells confirmed the three patient groups and individual assignment (Fig.1b). PRRX1-negative tumors were predominantly classified as Luminal A or B (65%), with a smaller fraction corresponding to HER2 (12%) or Basal (23%). This is compatible with previous observations showing that Prrx1-deficient PyMT tumors are highly differentiated19. In contrast, tumors with PRRX1-intermediate levels were more evenly distributed among Luminal B (22%), HER2 (33%), and Basal (45%), and all PRRX1-high tumors belonged to the Basal subtype.

To better understand the pattern of PRRX1 expression in human breast cancer, we analysed the METABRIC database5,24,25. Interestingly, almost 20% of patients show amplification of the PRRX1 locus (Supplementary Fig. 1a). Analysis of the Integrative Clusters described in the original METABRIC manuscripts revealed a clear enrichment of IntClust8 in patients bearing amplifications (Supplementary Fig. 1b), which is expected because this subgroup bears amplification of 1q chromosomal regions where PRRX1 is located. However, IntClust4 showed even higher PRRX1 expression than IntClust8 (Supplementary Fig. 1c) and is considered to have a good prognosis24. Regarding molecular subtypes, Claudin-low tumors exhibited the highest PRRX1 expression (Supplementary Fig. 1d), consistent with the strong EMT activation observed in this subtype26,27. Most Claudin-low tumors fall within IntClust4 (Supplementary Fig. 1e), the cluster with the highest PRRX1 levels (Supplementary Fig. 1f).

Finally, we observed a non-linear correlation between PRRX1 levels and metastatic incidence. The PRRX1-intermediate group exhibited the highest metastatic burden, with around 78% of patients developing metastases, compared to 31% in the PRRX1-negative group and only 12% in the PRRX1-high group (Fig. 1d and Supplementary Table 4). These findings are aligned with previous observations in patients showing that low Prrx1 expression is associated with worse survival20,21. These results suggest that PRRX1 levels, rather than its mere presence or absence, may be a critical determinant of the metastatic potential of cancer cells.

In vivo genetic modeling of different Prrx1 levels in cancer cells reveals a non-linear correlation with metastatic burden

We generated an in vivo mouse model to study the impact of different Prrx1 levels in cancer cells, by using a combination of (i) MMTV-PyMT28 transgenic mice, which develop mammary carcinomas that progress to invasive and metastatic stages resembling human invasive breast cancer29, (ii) a previously generated Prrx1 cKO mouse line carrying either one (Prrx1 +/f) or two (Prrx1 f/f) floxed Prrx1 alleles 19, and (iii) the Krt14-Cre tdTomato reporter line (Supplementary Fig. 2a-b). The Krt14-Cre line tags mammary gland progenitor cells from early embryonic stages30 (Supplementary Fig. 2b). Thus, in the adult mammary gland, all epithelial cells (luminal and basal) are Tomato⁺, regardless of Krt14 expression (Supplementary Fig. 2b,c). In PyMT tumors, virtually all cancer cells (PyMT+) are tdTomato+ (≈99% in previous studies19), irrespective of Krt14 expression (Supplementary Fig. 2b,d,e), while stromal cells remain negative (Supplementary Fig. 2d). Consistent with our recent finding that luminal cells acquire basal-like traits during PyMT tumor progression19, almost 99% of K14⁺ cells are PyMT⁺, ruling out non-transformed basal-like cells, as the oncogene is active only in the luminal compartment (Supplementary Fig. 2b,e). Altogether, these data confirm that our model is highly sensitive and specific, enabling Prrx1 deletion in all cancer cells before tumor onset without affecting other populations.

Aware that Prrx1 has different isoforms31, our conditional knockout (cKO) design targets exon 2, which is shared by both isoforms and encodes the critical homeobox domain. We previously showed that this strategy effectively abolishes the expression of both isoforms in embryonic fibroblasts19. Furthermore, heterozygous cells show similarly reduced levels of both isoforms19. To confirm that this is the case when combined with PyMT mice and to assess Prrx1 levels in the different allelic combinations in vivo, we analysed both total Prrx1 expression and that of the two isoforms (Long and Short) in FACS-sorted cancer cells from our experimental groups. We confirmed that Prrx1 +/f cancer cells expressed approximately half levels of Prrx1 compared to Prrx1 +/+ cells and that both isoforms behave similarly (Fig. 2a), ruling out the possibility of differential isoform effects in our model. In addition, we found that Prrx1 levels did not affect primary tumor burden (Fig. 2b) or the overall cell proliferation rate (Fig. 2c).

Consistent with our observations in breast cancer patients, mice bearing Prrx1 +/f tumors exhibited the highest metastatic burden, driven by an increase in the number and the size of metastatic foci. In contrast, Prrx1 f/f mice showed a significant reduction in metastatic burden compared to controls (Fig. 2d) in agreement with our previous findings19. Given that Prrx1 promotes cancer cell invasion19, we also assessed tumor invasiveness. In the PyMT model, luminal tumor cells progressively acquire basal-like traits, including de novo Krt14 expression, coinciding with activation of an invasive EMT program19, and Krt14 reliably marks the invasive and disseminative population32. Therefore, we combined peripheral enrichment of Krt14 expression with the presence of disorganized laminin borders, a well-established morphological hallmark of invasive behavior, to accurately delineate invasive areas. Compatible with our data on metastatic burden, Prrx1 +/f tumors displayed a significant increase in invasive areas with respect to control tumors, while Prrx1 f/f tumors showed a marked reduction (Fig. 2e). This is in line with our previous data showing that Prrx1 is associated with the invasive EMT phenotype19, and consistently, we detect Prrx1/K14 double positive cancer cells in these areas (Fig. 2f). Our data from Prrx1 genetic breast cancer mouse models reveal a non-linear relationship between Prrx1 levels and metastasis. Together with patient data, our findings in mouse models indicate that while a minimum threshold of Prrx1 is required for invasion and metastasis, high levels hinder metastatic progression.

Spatial transcriptomics identifies a molecular signature of cancer cells in the invasive areas

To better characterize the invasive population, we performed single-cell spatial transcriptomics on paraffin-embedded sections from Prrx1 +/+ primary tumors. With a custom panel of 264 genes designed to identify major cell types in PyMT tumors and capture EMT heterogeneity in cancer cells (Fig. 3a), we identified 1.794 million high-quality cells (Supplementary Fig. 3a,b), which were grouped into 18 transcriptionally distinct clusters with consistent representation across samples (Supplementary Fig. 3c). Using supervised annotation based on previously generated single-cell PyMT transcriptomics data19, we classified cells into five major types, encompassing cancer cells and associated populations, each characterized by specific marker expression (Fig. 3b,c). As expected, cancer cells were arranged in dense islands, surrounded by tumor microenvironmental populations, some of which infiltrated the tumor (Supplementary Fig. 3d).

To analyze the invasive regions, we manually identified areas of high keratin 14 (Krt14) expression at tumor margins (Fig. 3d), applying a strategy similar to the quantification of invasive areas (Fig. 2e). Comparing the expression profiles of cancer cells at invasive edges (Fig. 3e) with those located elsewhere revealed that invasive cells were significantly enriched in invasion-related genes (Tnc, Mmp2, Mmp3) and exhibited a basal-like expression program (Krt5, Trp63) (Fig. 3f and Supplementary Fig. 3e), while showing lower levels of injury response and interferon-related gene transcripts (Fos, Klf4, Mx1, Ifit1) (Fig. 3f). Invasive cells exhibited an EMT transcriptional profile (Fig. 3g), associated with transcription factors such as Snail2, a key regulator of the basal program33, Fibrillin-1 (Fbn1) and, importantly, Prrx1 (Fig. 3h). Visual inspection of the spatial transcriptomics data confirmed that this gene program is activated in cancer cells. As several of the detected genes, including Col1a1, Postn and Tnc, are typically highly expressed in stromal populations, we carefully verified that these results were not biased by transcript misassignment or imperfect cell segmentation at tumor borders (Supplementary Fig. 4). Overall, our spatial transcriptomics analysis revealed a molecular signature of cancer cells located at the invasive areas, highlighting a program compatible with that described along the Prrx1-dependent EMT invasive trajectory in scRNA-seq datasets19.

Intermediate levels of Prrx1 promote an invasive, highly proliferative metastatic state

To further dissect how different Prrx1 levels shape the transcriptional landscape of invasive cells, we performed scRNA-seq (10X Genomics) from 3 Prrx1 +/+, 3 Prrx1 +/f, and 4 Prrx1 f/f tumors (see methods). Unsupervised clustering identified nine major cell types and cancer cell states, each defined by bona fide lineage markers (Supplementary Fig. 5a-d). Focusing on cancer cells, we subset them and identified 28 subpopulations (Fig. 4a).

To identify the invasive population, we developed a Spatial Transcriptomics-derived Invasive Signature (STINS) using the differentially regulated genes (log2FC > 0.25 and Pvalue < 0.05) in the cancer cells residing in the invasive areas compared to those located elsewhere (Fig. 3f; See methods) and applied it to the scRNA-seq data, identifying clusters 20 and 25 as enriched for invasive features (Fig. 4b). These invasive cell populations present the highest levels of EMT activation in the tumor (Fig. 4c), with downregulation of the epithelial signature and upregulation of mesenchymal genes in cluster 25, exemplified by Epcam and Vim, respectively (Supplementary Fig. 6a). This is further supported by protein expression (Supplementary Fig. 6b). Clusters 20 and 25, are enriched in a recently described pro-invasive gene signature (PINGs)19 able to discriminate between invasive and non-invasive human breast cancer cell lines (Fig. 4d). Interestingly, these populations closely align with the final stages of the recently described invasive embryonic-like EMT trajectory in PyMT tumors (Supplementary Fig.7a,b). Notably, Prrx1 expression was only detected in cluster 25 (Fig. 4e).

Although we did not see differences in the distribution of the clusters (Supplementary Fig. 6c), we observed significant differences in the transcriptomic profile of cells in cluster 25 between the different genotypes. Comparative analysis of cluster 25 revealed that while the complete loss of Prrx1 (Prrx1 f/f) markedly reduced mesenchymal features and invasive signatures, Prrx1 +/f cells retained an EMT phenotype similar to the Prrx1 +/+ cells despite their lower level of Prrx1 expression (Fig. 4f). In addition, these findings indicate that intermediate Prrx1 levels are sufficient to sustain an EMT-driven invasive program, compatible with a partial and invasive EMT phenotype, previously described as highly metastatic15,16,18.

To better understand the increased metastatic potential of Prrx1 +/f tumors, we compared the expression profiles of cells in cluster 25 in Prrx1 +/+ and Prrx1 +/f conditions. Pathway enrichment analysis revealed that genes upregulated in Prrx1 +/f cells were primarily associated with proliferation and cell cycle regulation (Fig. 4g). This was further confirmed using cell cycle scoring (Reactome), which indicated that Prrx1 +/f invasive cells are significantly more proliferative than their Prrx1 +/+ counterparts (Fig. 4h). Immunofluorescence quantification of cycling cancer cells at the invasive edge in Prrx1 +/+ and Prrx1 +/f tumors validated this observation (Fig. 4i), compatible with the enlarged invasive areas observed in Prrx1 +/f tumors (Fig. 2e). Importantly, this effect on proliferation is restricted to the invasive cell population. As expected from an effect dependent on Prrx1 dose, we did not observe significant changes in primary tumor size and global proliferation (Fig. 2b,c), as over 95% of tumor cells do not express Prrx1 (Fig. 4e)

To identify potential regulators of changes in proliferation in cluster 25, we examined the expression of cell cycle-related genes in invasive cells and found a notable enrichment of Cyclin D2 and CDKN2 family members when comparing cells in cluster 25 with the rest (Fig. 4j). Given that CDKN2 proteins are key inhibitors of the cell cycle34, and cell cycle arrest is one of the main characteristics of dormancy, we investigated whether the Prrx1-positive invasive cell population also expressed dormancy-associated genes. Using a breast cancer dormancy signature recently derived by comparing cycling and non-cycling metastatic PyMT cells, able to differentiate between the dormant and non-dormant version of the same cells (D2.A1 vs D2.OR cell lines35), we found that our invasive clusters had the highest dormancy scores in the tumor, cluster 25 in particular (Fig. 4k). Immunofluorescence analysis of Mme expression, a marker of the dormancy program, confirmed high expression in invasive areas (Fig. 4k), indicating that dormancy and invasion are coregulated in tumors. Notably, Prrx1 +/f cells in cluster 25 exhibited a significantly lower dormancy score than their Prrx1 +/+ counterparts (Fig. 4k), suggesting that intermediate Prrx1 levels allow invasive cells to avoid the activation of dormancy. To assess whether these differences in proliferation and dormancy at the invasive edge could be translated to the metastatic site, we quantified the proportion of proliferating cancer cells in lung metastases, a common approach to evaluate dormancy at the metastatic site36,37. Metastases derived from Prrx1 +/f tumors displayed a significantly higher fraction of cycling cells compared to those from Prrx1 +/+ tumors (Fig. 4l), consistent with the previously observed increased size of Prrx1 +/f metastatic foci (Fig. 2d). Finally, we aimed to determine whether this state was also associated with enhanced tumor-initiating capacity using CD90, a well-characterized marker of metastatic-initiating cells in the PyMT model38. FACS analysis confirmed that the Epcam⁻ cancer cell population - and cluster 25 shows the lowest Epcam expression in the tumor (Supplementary Fig. 6a) was enriched in CD90⁺ cells in Prrx1 +/f tumors compared to the same population in control Prrx1 +/+ tumors, suggesting that intermediate Prrx1 levels promote metastatic-initiating capabilities (Fig. 4m). Interestingly, this finding aligns with previous observations in cultured human cell lines, where PRRX1 downregulation enhanced sphere formation increasing the proportion of CD44⁺/CD24⁻ cells20. Overall, these findings suggest that while high Prrx1 expression maintains invasive cells in a dormant, low-proliferative state, intermediate levels sustain invasive potential while allowing proliferation, facilitating metastatic colonization. This mechanism could explain the higher metastatic incidence observed in patients bearing breast tumors with intermediate PRRX1 levels.

Chromatin accessibility analysis establishes Prrx1 as a direct regulator of proliferation and dormancy in the invasive population

To further investigate how Prrx1 could regulate cell cycle progression and dormancy, we profiled chromatin accessibility after snATAC-seq and analyzed chromatin dynamics in Prrx1 control tumors (Fig. 5 and Supplementary Figs. 8 and 9). Unsupervised clustering of open chromatin data identified six main clusters (Supplementary Fig. 9a), with a consistent cell distribution across samples (Supplementary Fig. 9b). We used gene activity of bona fide cell type markers to annotate clusters into major cell types (Supplementary Fig. 9c-d). As for the scRNA-seq analysis, we focused on the cancer cells and identified 16 distinct clusters (Fig. 5a). To identify invasive cells, we imputed gene expression in sATAC-seq cells using chromatin accessibility data and scRNA-seq as a reference (see Methods) and applied STINS. As observed in the scRNA-seq data, two clusters, 14 and 15, emerged as the invasive populations (Fig. 5b), with cluster 15 exhibiting higher Prrx1 imputed expression (Fig. 5c), mirroring cluster 25 in the scRNA-seq analysisAs observed in the scRNA-seq dataset, we also observed a correlation between chromatin states and our recently described EMT trajectories, indicating that distinct transcriptional EMT states are supported by specific chromatin profiles (Supplementary Fig. 10a-b).

Given the nature of Prrx1 as a transcription factor, we hypothesized that it could be a master regulator of the invasive and metastatic cells. Transcription factor binding motif analysis showed that Prrx1 binding sites are more accessible in cluster 15 compared to the rest of the clusters (Fig. 5d). Further analysis revealed that several accessible regions in cluster 15, located near the promoters of genes encoding cell-cycle regulators (cyclins and cyclin inhibitors) and dormancy-related genes³⁴ (Gas6, Mme), respectively enriched in the invasive cell population in the scRNA-seq dataset (Fig. 4j) and part of the dormancy signature (Fig. 4k), are positive for Prrx1 binding motifs as assessed with the 'Find Individual Motif Occurrences' tool (FIMO; Fig. 5e,f).

These results indicate that Prrx1 is not only specifically expressed in the invasive population but also that it plays an important regulatory role, as it can directly bind promoters of genes associated with cell-cycle control and dormancy, providing mechanistic evidence for the differences in metastatic potential observed at different Prrx1 levels.

Engineered cell line models reveal dose-dependent Prrx1 regulation of invasion, proliferation and dormancy

To examine whether this regulatory pattern was conserved in human, we analyzed whole-genome transcriptomes of human breast cancer cell lines39 ordered and grouped by EMT status19 (Fig. 6a). Firstly, we compared PRRX1 expression across EMT groups and confirmed that the mesenchymal group shows the highest levels (Supplementary Fig. 11a), consistent with PRRX1 being enriched in mesenchymal cells and expressed in PyMT tumors in cells with the highest EMT score19. Hybrid E/M groups displayed intermediate expression when compared to epithelial and mesenchymal states (Supplementary Fig. 11a).

When scoring the different gene signatures in this dataset, a clear correlation emerged between PRRX1 expression, mesenchymal character, invasion (PING) and dormancy signatures (Fig. 6a). To assess its predicted regulatory role, we selected BT549 cells, as they express high levels of PRRX1 and of dormancy and cell cycle inhibitor genes (Fig. 6a), the closest to cells in cluster 25 from the scRNA-seq data in PyMT tumors . In addition, we had previously shown that PRRX1 downregulation is necessary for BT549 cells to metastasize while it reduced the invasive capacity of the BT-549 cells20, making them an ideal model for assessing the impact of different PRRX1 downregulation on cell cycle regulators and dormancy-related genes. Downregulation of PRRX1 by RNA interference (BT549 shPRRX1 cells) led to an increase in the expression of the CCND1 cyclin and a reduction in that of cyclin inhibitors of the CDKN2 family and dormancy-associated genes (GAS6, MME, OGN) (Fig. 6b). In addition, analysis of published data on PRRX1 downregulation in the HS578T human breast cancer cell line40 revealed a clear upregulation of CCND2, another Cyclin D family member, while changes in CCND1 were less evident, along with downregulation of MME, a dormancy-associated gene (Supplementary Fig. 11b). Finally, the analysis of Chip-seq dataset in three different human breast cancer cell lines with different PRRX1 levels40, including HS578T, shows that, at least one of them, PRRX1 can directly bind the promoters of cell cycle regulators, dormancy and invasion genes in human breast cancer cell lines (Fig. 6c and Supplementary 11c). These findings confirm that PRRX1 is an activator of cell cycle inhibitors and dormancy-associated genes in human breast cancer and provide an explanation for the requirement of PRRX1 downregulation for metastatic outgrowth.

To mechanistically model graded Prrx1 expression we engineered cell-line variants expressing low (CTR), intermediate, or high levels of Prrx1 (Fig. 6d). The difference between intermediate and high expression was approximately twofold, closely matching the allelic differences observed in our mouse model. We selected MDCK-NBL2 epithelial cells as the parental line because of their pronounced phenotypic plasticity, well-characterized EMT response to TGF-β19,41–44, and their ability to recapitulate the sequential activation of EMT transcription factors observed in neural crest cells and in the PyMT model, including the late induction of Prrx119. Importantly, the MDCK system allows us to study the effects of increasing Prrx1 dosage on cell-state transitions and to directly assess its relationship with invasion and proliferation. Both intermediate- and high-Prrx1-expressing clones displayed similar and significant invasiveness compared to CTR cells (Fig. 6e). In contrast, proliferation assays revealed a marked reduction in growth only in Prrx1-high clones (Fig. 6f). Gene expression analysis mirrored these functional outcomes. Prrx1-high clones exhibited reduced expression of CCND1 and increased levels of CDKN2B/C, whereas intermediate clones remained closer to CTR cells (Fig. 6g). Consistent with this, dormancy-associated genes (GAS6, MME) were strongly upregulated in Prrx1-high clones (~100-fold) but downregulated in intermediate clones (Fig. 6g). Collectively, these data indicate that high Prrx1 levels promote invasion while suppressing proliferation and inducing dormancy, whereas intermediate Prrx1 expression enhances invasion with minimal impact on proliferation, thereby maximizing metastatic potential. These findings recapitulate our in vivo observations and support a dose-dependent, hormetic model of Prrx1 function regarding the regulation of metastatic potential.

Patient stratification by combining invasion and proliferation signatures reveals a high-risk breast cancer subgroup with intermediate Prrx1 levels

We next wanted to assess whether the relationship between PRRX1 levels and invasion, proliferation, and dormancy observed in the mouse tumors and human cell lines was conserved in human breast cancer. We analyzed the human breast cancer METABRIC database24,25 and found that PRRX1 expression positively correlates with signatures associated with invasion (PINGs) and dormancy, and negatively correlated with proliferation (Fig. 7a). Notably, this relationship extends to genes identified as putative PRRX1 targets in our previous analyses, including MMP2, ADAMTS2, FBN1, COL1A1, CCND1, CDKN2B/C, CFH, GAS6, OGN, MME (Supplementary Fig.12).

To evaluate whether the interplay between invasion and proliferation has prognostic value, we stratified METABRIC patients5,24,25 into four groups (see Methods) based on their scores for the invasion (PING)19 and proliferation (Cell_Cycle_Mitotic_R-HSA-69278) signatures, both of which were found to be dependent on PRRX1 levels in the invasive population in our genetic mouse models (Fig. 4f and 4h).

Interestingly, the impact of high invasion on survival was contingent on proliferation status, patients with highly invasive but low-proliferative tumors showed better prognosis than the rest of the groups (Fig. 7b and Supplementary Table 4). More importantly, this classification revealed that patients with tumors scoring high for both invasion and proliferation exhibit significantly worse overall survival (OS), particularly during the first months after diagnosis (Fig. 7b and Supplementary Table 4), compatible with our results modeling different Prrx1 levels in the mouse tumors and with our proposal that metastatic competence arises from the coexistence of both programs, invasion and proliferation. Importantly, even when stratified by treatment, a very similar pattern emerged across all conditions (Supplementary Fig. 13).

Given the established link between invasion and proliferation, we next investigated whether our patient stratification reflected differences in dormancy scores and found that, as expected, patients from the group with the poorest survival, highly invasive and highly proliferative tumors showed a markedly lower dormancy score (Fig. 7c). Interestingly, this group of aggressive tumors show lower PRRX1 expression levels than the invasive, low-proliferative group (Fig. 7d). These data fit with the non-linear metastatic burden observed in TMAs from patients with different PRRX1 expression levels (Fig. 1). Taken together, these findings indicate that human breast cancer patients with tumors showing intermediate levels of PRRX1 in cancer cells exhibit a unique combination of invasive and proliferative traits, which defines a particularly aggressive subset of tumors. Moreover, our stratification approach could serve as a valuable tool to identify high-risk patients who may benefit from therapeutic strategies targeting both invasion and proliferation.

Methods
Breast Cancer Tissue Microarray (TMA) Analysis
A total of 94 infiltrating grade 3 ductal female breast carcinoma (IDC) samples obtained from the archives of the MD Anderson Cancer Center, Madrid, Spain (MD Anderson Foundation Biobank, record number B.0000745, ISCIII National Biobank). Patients underwent surgery between 2010 and 2011, with a mean age at surgery of 57.07 years (range: 44–81 years). According to the TNM Classification staging, 32 tumors were stage I, 31 stage II, and 31 stage III–IV. Histological and immunohistochemical studies were performed on formalin-fixed, paraffin-embedded tissue samples57. Two 1-mm cylinders were taken from each tumor block to generate TMA blocks58. Inter-observer concordance between pathologists was high, exceeding 90%, in line with standard practice for immunohistochemical evaluation in breast pathology. All cases were independently assessed by two experienced pathologists. In instances where substantial discordance was observed, the immunostaining was repeated and evaluated on full tissue sections rather than on TMA cores, allowing for comprehensive reassessment and consensus scoring. The study followed standard ethical procedures outlined in Spanish regulations (Ley de Investigación Orgánica Biomédica, July 14, 2007, and Real Decreto Biobancos 1716/2011) and was approved by the ethics committee of the MD Anderson Cancer Center Madrid, Spain.

Animal Experiments
Mice were fed ad libitum. Housing and experimental procedures were conducted in strict compliance with the European Community Council Directive (89/609/EEC) and Spanish legislation. Ethical protocols were approved by the CSIC Ethical Committee and the Animal Welfare Committee of the Institute of Neurosciences. Animals for experiments were selected by genotype, and no randomization or blinding was performed. Animals were housed under SPF conditions at the RMG Animal House (ES-03-119-0001001 SEARMG). Mouse experiments were performed in the MMTV -PYMT model crossed with a Rosa-LSL-tdTomato reporter line, purchased from JAX Mice (The Jackson Laboratory), expressing tdTomato upon Cre-mediated recombination. Cre recombinase is expressed under the control of the KRT14 promoter (Tg(KRT14-cre)1Amc/J; 004782). Mice were backcrossed in the FVB background for at least ten generations (99.9% FVB) The generated line was further crossed with Prrx1 conditional mutant mice (Prrx1em1An)19 to give rise to heterozygous (Tg(MMTV-PyMT)634Mul; Tg(KRT14-cre)1Amc; Prrx1em1An /+ named as Prrx1 +/f) and homozygous (Tg(MMTV-PyMT)634Mul; Tg(KRT14-cre)1Amc; Prrx1em1An named as Prrx1 f/f ) conditional Prrx1 mutants. Considering that the study focuses on breast cancer, only female PyMT mice were used. Mammary gland carcinomas were collected from 14 to 16-week-old mice, upon reaching the maximum tumor size or burden permitted by ethical guidelines. Maximal tumor size/burden was defined as equal to or greater than 1,500 mm3 per tumor or equal to or greater than 3,000 mm3 total tumor burden following the ethical protocols approved by the CSIC Ethical Committee and the Animal Welfare Committee of the Institute of Neurosciences (protocols 2015/VSC/PEA/00211 and 2019/VSC/PEA/0218). Maximum size was never exceeded. Health state, tumor size and burden were monitored weekly.

Immunofluorescence (IF)
Tumor and Lung Sections: Sample preparation and imaging: 3–5 µm paraffin-embedded sections were dewaxed, and protein epitopes unmasked by immersion in 95°C preheated citrate buffer (pH 6.0) or Tris-EDTA buffer (pH 9.0) for 20 minutes. O.C.T.™ sections (5 µm) were washed three times in PBS for 5 minutes. For IF staining, slides were placed in a humidified chamber and blocked/permeabilized for 1 hour with IF blocking buffer (5% normal goat or donkey serum, 1% BSA, and 0.2% Triton X-100 in PBS). Blocking solution was replaced with the primary antibody diluted in cold IFBB, and the chamber was incubated overnight at 4°C. Slides were washed 3× for 10 minutes in PBS, incubated for 1 hour with secondary antibodies and DAPI, then washed again 3× for 10 minutes in PBS.Slides were mounted using anti-fade mounting medium (DAKO). Primary and secondary antibodies used are shown in Supplementary Table 1. Images were captured using a Leica DMR microscope for brightfield images and an Olympus FV1200 microscope for confocal imaging.

Whole Lungs: TdTomato+ metastases were visualized by IF on whole lungs and cleared following the iDISCO+ protocol59. Images were acquired using an UltraMicroscope II (LaVision BioTec), analyzed, and reconstructed in 3D using Vision4D (Arivis) Image Analysis Software. Metastatic burden was assessed using Imaris software (version 9.3.1; BitPlane). The “Surface” function was used for segmentation (tdTomato signal), and volumetric data were extracted. To avoid false positives, a detection cutoff of 90,000 µm³ was identified as the minimum volume for high-confidence object identification.

Mammary Tumor and Lung Sample Collection and Preparation
Tumor and whole lung samples were fixed in 4% PFA overnight at room temperature (RT) for embedding in paraffin or for 4 hours at 4°C for embedding in O.C.T.™ (Sakura). The next day, samples were washed several times and stored in PBS + Azide overnight at 4°C. Pre-fixed samples were embedded in paraffin or O.C.T.™ for further sectioning and collection on SuperFrost Plus microscope slides.

Cancer Cell FACS Sorting and Analysis
Mammary gland carcinomas were collected from 14 to 16-week-old female mice. Harvested tissues were manually minced using sterile scalpels and finely chopped with a McIlwain Tissue Chopper (TED PELLA, INC). Minced samples were incubated in 2.5 ml of digestion buffer (PBS supplemented with 1 Wünsch units of TH Liberase/ml and 25 μg/ml DNase I (Roche) at 37°C for 75 minutes in an Incubator Microplate Shaker (VWR). Tumor fragments were pipetted up and down every 5 minutes to achieve tissue dissociation into single-cell suspensions, and disaggregation was monitored under a microscope. Disaggregated samples were neutralized with 15 mL of breast medium supplemented with 25 μg/ml DNase I and passed through 70 μm and 40 μm filters (BD Falcon). Cells were centrifuged, and pellets were resuspended in 5 mL of red blood cell lysis ACK buffer (Ammonium-Chloride Potassium) for 5 minutes at room temperature (RT) with gentle rotation. After red blood cell removal, cancer cells were neutralized, resuspended in 0.5 ml FACS buffer with DAPI per 1M cells, and directly used for FACS to sort Tomato⁺ cells. Cell sorting was performed using a BD FACSAria III flow cytometer, and 300,000 cancer cells were sorted. The sorted cells were centrifuged at 5,000 rpm for 3 minutes and lysed using the lysis buffer from the illustra RNAspin Mini Kit (GE Healthcare) following the manufacturer’s instructions.
For immunolabeling, after neutralization of the red blood cell lysis buffer, cells were resuspended in 200 µl of FACS buffer per 1M cells. Fc receptor blocking was performed by adding Palex-BioLegend 156604 (1:50) and incubating for 15 minutes under agitation at 4°C. Fluorescence-labeled antibodies were then added, and cells further incubated for 30 minutes under agitation at 4°C. Finally, cells were resuspended in 200 µl of FACS buffer with DAPI per 1M cells and analyzed using the BD FACSAria III flow cytometer and BD FACSDiva (BD Bioscience) software.

Total RNA Extraction, cDNA Synthesis, and RT-qPCR
For gene expression analysis, RNA was extracted using the Illustra RNAspin Mini Kit (GE Healthcare). Reverse transcription was performed using the Maxima First-Strand cDNA Synthesis Kit (ThermoFisher Scientific). RT-qPCR was carried out using Fast SYBR Green Mastermix (Applied Biosystems) on a Step One Plus machine (Applied Biosystems), following the manufacturers’ instructions. Relative RNA expression levels (relative fold change) were calculated using the 2^−ΔΔCt formula. Quantitative RT-qPCR primers were ordered to IDT and sequences are listed in Supplementary Table 2.

Patient clustering based on PRRX1 expression
To quantify PRRX1 expression of patients, we computed two metrics for each tumor sample, the percentage and the mean intensity of PRRX1 expression in cancer cells. These values were extracted from single-cell segmentation of IF data using CellProfiler and CellAnalyst softwares60,61. The resulting two-dimensional feature space was standardized using StandardScaler from the sklearn.preprocessing module. Unsupervised clustering was then performed using the K-means algorithm (KMeans from sklearn.cluster). The resulting cluster labels were assigned to each patient and used to color-code data points in the scatter plot.

Spatial transcriptomics analysis
Spatial Transcriptomics Panel Design: We designed a MERSCOPE panel of 265 genes (259 genes in the main panel and 6 in the sequential panel) using the Vizgen Gene Panel Design Portal. The genes included are shown in Supplementary Table 3. The gene selection was based on two criteria: Firstly, genes that cover the main cell types and cell states (primarily cancer cell states) identified in previously published scRNA-seq data of PyMT tumors19. Secondly, genes involved in general signaling pathways such as Wnt, Hedgehog and Hippo pathways, plus cytokine-cytokine receptor interactions, prioritizing those detected in the scRNA-seq data mentioned above. To avoid optical crowding during transcript decoding, we analyzed public bulk RNA-sequencing data from PyMT tumors and candidate genes with very high expression levels were excluded. For each gene, a panel of 30 encoding probes was designed by Vizgen using a proprietary algorithm. The final gene panel also included 30 blank probes serving as negative controls.

Spatial Transcriptomics Sample Preparation and Data Acquisition: Formalin-fixed paraffin-embedded (FFPE) tissue from 3 Prrx1 +/+ PyMT primary tumor was sectioned to a thickness of 5 μm following standard histology protocols. Sections were floated for 5 minutes in a 45°C water bath, mounted on fiducial bead-coated FFPE MERSCOPE slides (Vizgen, #20400100), dried at room temperature in a fume hood for 1 hour, baked at 55°C for 15 minutes, dried at room temperature for 2 hours, and stored at -20°C overnight. Samples were prepared following the Vizgen MERSCOPE® User Guide for FFPE Tissue Sample Preparation, beginning at the cell boundary staining step. To remove lipids and proteins that interfere with imaging, 5 ml of Clearing Premix (Vizgen #20300003) were mixed with 100 μl of Proteinase K for each sample. The samples were incubated at 47°C in a humidified incubator overnight. Samples in Clearing Solution were photobleached for 3 hours to remove background fluorescence using the MERSCOPE Photobleacher (Vizgen #1010003). Following photobleaching, the samples were incubated overnight at 37°C, with the clearing solution replaced with 200 μL of Clearing Premix containing 50 μL of Proteinase K for an additional 4 hours at 37°C. Finally, the samples were incubated for 48 hours with the customized gene panel described above. Samples were stored in Clearing Solution at 37°C in the incubator prior to imaging for up to one week. Imaging was performed on the MERSCOPE platform according to the MERSCOPE Instrument User Guide. Each field of view was imaged in seven 1.5 μm-thick z planes at 60× magnification. Images were processed and decoded to identify RNA spots with xyz coordinates and gene assignments using Vizgen's Merlin software.

Cell Segmentation: To improve the MERSCOPE onboard cell segmentation, we re-segmented the data using Cellpose2 v2.2.362 implemented in the vizgen post-porcessing tool (VPT) v1.2.2 framework (https://github.com/Vizgen/vizgen-postprocessing). Cyto2 model provided with the Cellpose2 was used for the segmentation. The z_layers was set to 3 and the three image channels i.e. Cellbound2, Cellbound3, and DAPI were mapped to the red, green, and blue channels, respectively. In the preprocessing, gaussian blur (kernel size = 11) was used to smooth the image boundaries. Intensity was normalized using contrast-limited adaptive histogram equalization (CLAHE) with a clip limit of 0.01 and a filter size of 100 × 100 pixels. For the segmentation, DAPI was used as a nuclear channel and object diameter threshold was set to 70 pixels to specify approximate cell size. The diameter threshold was derived from the manual inspection of random areas of MERSCOPE onboard segmentation in cellpose2 visualizer. A flow threshold of 0.95 and a cell probability threshold of -5.5 were applied, with a minimum mask size of 500 pixels to filter the noise. Post-segmentation, polygon smoothing, and cell boundaries simplification were performed with a tolerance of 2 pixels and a smoothing radius of 10 pixels, retaining only polygons with a minimum area of 500 pixels. To harmonize the results, entities were fused with a minimum distance of 1 pixel between them. Finally, partition-transcripts utility from VPT framework was used to obtain the cell by gene matrix using re-segmented data and derive-entity-metadata was used to derive the cell metadata which includes, the cell coordinates, volume, and transcript counts for each cell. The generate-segmentation-metrics was used to generate the report which manually inspected to ensure the quality of re-segmentation of each sample.

Quality control, filtering, and data integration: The cell-by-gene matrix, along with cell metadata, was imported into the Python environment and converted into an AnnData object using Scanpy v1.9.363. Thresholds for gene expression and cell volume were determined by manually inspecting their distributions across each sample. Cells with fewer than 15 gene expression counts, volumes smaller than 50 µm3, or larger than 5000 µm3 were removed from downstream analysis which resulted in 1.794 million total number of cells. AnnData objects from three different samples were merged, and library normalization was performed using the scanpy.pp.normalize_total function. Then the count matrix was log-transformed using the scanpy.pp.log1p function and scaled using the scanpy.pp.scale function (max_value=10). Principal component analysis (PCA) was performed on the scaled data using scanpy.tl.pca function with default parameters. Finally, sample-level data integration was carried out using Harmony64, through scanpy.external.pp.harmony_integrate function.

Dimensionality reduction, clustering, and cell type annotation: The neighborhood graph was constructed using the top 20 Harmony components with the scanpy.pp.neighbors function (n_neighbors=10). The data were represented in low dimensional space using Uniform Manifold Approximation and Projection (UMAP) using scanpy.tl.umap function (min_dist=0.5 and spread=1.0). The initial round of Leiden clustering was performed using the scanpy.tl.leiden function with a resolution of 0.5. After removing low-quality clusters, a second round of clustering was conducted at a resolution of 1.2, following the same pipeline, starting with re-normalizing the count matrix which resulted in 18 total number of clusters (1.793 million cells). Differential gene expression analysis was carried out using the Wilcoxon test (scanpy.tl.rank_genes_groups), and the expression of bona fide marker genes was used to annotate the cell clusters and expression were represented as dotplot to show their specificity using scanpy.pl.rank_genes_groups_dotplot function.

Identification of invasive edges MERSCOPE data: The cell metadata along with cell type annotation was extracted from the AnnData object and imported in MERSCOPE visualizer. The invasive areas were drawn manually by looking for peripheral tumor edges enriched in Krt14 expression. Cancer cells were categorized as part of invasive edge or bulk of the tumor. The metadata for cancer cell annotated as part of invasive edge or bulk of the tumor was exported and used for the downstream analysis.

Differential gene expression between cancer cells in invasive edge and bulk of the tumor: Differential gene expression analysis was performed to compare the cancer cells located at the invasive edge and the bulk of the tumor. The raw count matrix and associated metadata were imported into the R v.4.3.3 (https://www.r-project.org/) statistical framework and converted into a Seurat (v.5.1.0) object using the CreateSeuratObject function. The cancer cells were subsetted from the Seurat object using subset function and normalized to the library size using the NormalizeData function. The FindMarkers function was used to identify differentially expressed genes between cancer cells at the invasive edge and in the tumor bulk. The resulting statistics were visualized as a volcano plot using ggplot2. Genes upregulated in the invasive edge were identified as the invasive gene signature. The pathway enrichment analysis was performed for the differentially regulated genes (log2FC > 0.25 and Pvalue < 0.05) using EnrichR65 against MsigDB and ARCHS4 TFs database. The same gene set was used to create spatial transcriptomics-derived invasive signature (STING).

Single cell RNA-seq analysis
Single-Cell Preparation: Mammary gland carcinomas were collected from 14- to 16-week-old female mice. Harvested tissues were manually minced using sterile scalpels and finely chopped with a McIlwain Tissue Chopper (TED PELLA, INC). Minced samples were incubated in 2.5 mL of digestion buffer (PBS supplemented with 2.5 Wünsch units of TH Liberase/mL and 25 μg/ml DNase I (Roche) at 37°C for 45 minutes in an Incubator Microplate Shaker (VWR). Tumor fragments were pipetted up and down every 5 minutes to achieve tissue dissociation into single-cell suspensions, and disaggregation was monitored under a microscope. Disaggregated samples were neutralized with 15 ml of breast medium supplemented with 25 μg/ml DNase I and passed through 70 μm and 40 μm filters (BD Falcon). Cells were centrifuged, and pellets were resuspended in 5 mL of red blood cell lysis ACK buffer (Ammonium-Chloride Potassium) for 5 minutes at room temperature (RT) with gentle rotation. Following another centrifugation, pellets were resuspended in 5 mL of FACS buffer/5 mM EDTA and passed through a 20 μm filter. To assess single-cell preparation, 20 μl of cell suspension was mixed with an equal volume of Trypan Blue and incubated for 5 minutes to evaluate cell number, viability, and debris content. Cells were adjusted to a concentration of 1,000 cells/μl with cold FACS buffer and used directly for GEM (Gel Bead-In Emulsions) preparation.

Single-Cell GEM and cDNA Library Preparation: 3 Prrx1 +/+, 3 Prrx1 +/f and 4 Prrx1 f/f single-cell preparations were obtained from ten female mice, with each preparation comprising a mix of mammary gland carcinomas from the same mouse. Individual cell encapsulation for single-cell expression profiling and library preparation was performed using the Chromium Single Cell 3’ Library & Gel Bead Kit v3.1 (10x Genomics) and the Chromium Controller (10x Genomics). Samples were multiplexed using the Chromium Dual Index Kit TT Set A. Pooled libraries were loaded onto an Illumina Nextseq 2000 machine with a 100 bp read length.

Quality control, filtering, and data integration: scRNA-seq libraries for different genotypes were sequenced to recover 69436 total number of cells. The raw data quality was inspected using FastQC v0.11.9 (Babraham Institute). The reads were aligned to the mouse genome assembly GRCm38 (mm10; Ensembl reference annotation) and quantified using the 10X Genomics Cell Ranger v6.1.1 pipeline66 with default parameters. The CellRanger output was imported in the R v.4.3.3 statistical environment and converted to a Seurat (v.5.1.0) object using the CreateSeuratObject function. The low-quality cells were detected based on the number of genes (400-4000) and the percentage of mitochondrial genes (15%). The putative doublets were detected using scDblFinder v1.8.067. The low-quality cells and doublets were removed from the downstream analysis. The high-quality cells (n=50386) were then subjected to the library normalization using the NormalizeData function with default parameters followed by the high variable gene (HVGs) detection. In this analysis we integrated data from 4 similar stage PyMT tumors Prrx1 +/+ from19, which were prepared with 10X single cell 3' v2 chemistry. The top 4000 HVGs were selected for the Prrx1 f/f and Prrx1 +/f samples, whereas top 6000 HVGs were selected for Prrx1 +/+ samples. The detection of HVGs was performed separately for each sample. Then the intersection was taken for HVGs calculated for each sample per condition and finally a non-redundant set of genes was taken from the condition specific intersect and used for the sample integration. The gene expression of these HVGs was scaled using the ScaleData function and subjected to the PCA. Finally, the samples were integrated based on the top 30 principal components using Harmony64 implemented in the RunHarmony Seurat function.

Dimensionality reduction, clustering, and cell type annotation: A Shared Nearest Neighbor (SNN) Graph for the integrated Seurat object was built using FindNeighbors over the top 30 harmony components. The Uniform Manifold Approximation and Projection (UMAP), a nonlinear dimensional reduction was used for the data visualization calculated over the same top 30 harmony components using the RunUMAP function. The single cells were clustered using an SNN modularity optimization-based clustering algorithm implemented in FindClusters function (resolution: 0.035) which resulted in 9 clusters. First the gene counts for all the genes in RNA assay were Log-Normalized using NormalizeData function and the differential gene expression was calculated using FindAllMarkers with default parameters. The expression of bona fide markers for each cluster was represented as a dotplot using DotPlot function showing their specificity for a cluster.

Clustering of cancer cell subset: The cancer cell cluster was subsetted and re-clustered to better understand the transcriptional heterogeneity within cancer cell. Sample-specific HVGs were using the workflow described above in "Cell quality control, filtering, and data integration". In the initial round of clustering, we used the top 25 Harmony components with a resolution of 1.0. During this step, clusters with high expression of mitochondrial and ribosomal genes were identified and removed as poor-quality clusters. Subsequently, the cleaned cancer cells were re-integrated based on recalculated HVGs using Harmony. A second round of clustering was performed on the top 25 Harmony components, resulting in the identification of 28 distinct clusters at a resolution of 1.5. To annotate the identified cancer cell clusters, a gene signature scoring approach was applied. Briefly, the top 20 differentially expressed genes for each cancer cell cluster, were taken from19 and the gene signature score per cell was calculated using the AddModuleScore_UCell function from the UCell v2.0.0 package68 These scores, combined with known marker genes describing cancer cell states, were used to annotate the clusters. In similar fashion, the score was calculated for dormancy, inflammation, EMT, PINGs and spatial-derived invasive gene signature.

Condition specific differential gene expression analysis: To minimize technical variability introduced by different chemistries, for the condition comparison, we utilized only the Prrx1 f/f samples generated with the 10X single cell 3' v3 chemistry, as the Prrx1 +/f and Prrx1 f/f samples were produced using the same v3 chemistry. This sample composition was used while comparing the gene signatures score across different conditions in clusters of interest. To identify condition-specific differentially expressed genes within target clusters, a pseudobulk approach was employed using the muscat v1.16.0 package69. Genes with expression count greater than 5 in more than 3 cells were considered for the comparison. The gene counts for each sample within individual cell states were aggregated using the aggregateData function. Genes with a log2 fold change greater than 0.25 and a p-value less than 0.05 were considered as upregulated. These upregulated genes were then used for pathway enrichment analysis against the MsigDB database using EnrichR65.

Single cell ATACseq analysis
Single-Nuclei Preparation: Four mammary gland carcinomas from Prrx1 +/+ PyMT mice were collected from 14- to 16-week-old female mice and manually minced using sterile scalpels. Minced samples were transferred to a tissue homogenizer douncer prefilled with Solution D (20 mM Tris-HCl pH 7.5, 0.1% Tween 20, 0.25 M Sucrose, 25 mM KCl, 5 mM MgCl₂) and homogenized 8–10 times with rotation. The lysate was passed through 70 μm and 40 μm filters (BD Falcon) and transferred to a 15 mL tube. After centrifugation at 500 × g for 5 minutes at 4°C, the supernatant was removed, and the pellet was resuspended in 4 mL of Solution D. Two milliliters of Optiprep (STEMCELL) were added, and the sample was centrifuged at 1,500 × g for 10 minutes at 4°C. The pellet was resuspended in Wash Buffer (10 mM Tris-HCl pH 7.5, 10 mM NaCl, 3 mM MgCl₂, 2% BSA, 0.1% Tween 20, 1 mM DTT) and passed through a 40 μm filter. The sample was transferred to a 1.5 mL tube and centrifuged at 500 × g for 5 minutes at 4°C. Nuclei were resuspended in 200 μL of 0.1× lysis buffer (10 mM Tris-HCl pH 7.5, 10 mM NaCl, 3 mM MgCl₂, 2% BSA, 0.1% Tween 20, 1 mM DTT, 0.1% Nonidet P40, 0.01% Digitonin) and incubated on ice for 5 minutes. Following incubation, 1.3 mL of Wash Buffer was added, and the preparation was passed through a 20 μm filter. After centrifugation at 500 × g for 5 minutes at 4°C, nuclei were resuspended in 100 μL of 1× Nuclei Dilution Buffer (Chromium Next GEM Single Cell ATAC Kit v2). Ten microliters of the preparation were mixed with an equal volume of Trypan Blue, and nuclei number and viability were determined using the automated Countess II cell counter (Invitrogen).

Nucleus Transposition and Single-Nuclei GEM and DNA Library Preparation: Four single-nuclei preparations were obtained from four female mice, with each preparation comprising a mix of mammary gland carcinomas from the same mouse. Prepared nuclei suspensions were transposed following the guidelines of the Chromium Next GEM Single Cell ATAC Kit v2 (10x Genomics) and encapsulated using the 10x Chromium Controller. Samples were multiplexed using the Chromium Single Index Kit N. Pooled libraries were loaded onto an Illumina Nextseq 2000 machine with a 50 bp read length.

Quality control, filtering, and integration process: Libraries were sequenced to obtain 887.61 million reads in total. The quality of the sequenced reads was assessed using FastQC v0.11.9 (Babraham Institute). Sequenced samples were processed using the CellRanger-atac v2.1.0 pipeline (10x Genomics) and aligned to the GRCm38 (mm10) mouse reference genome (Ensembl annotation v99). MACS2 v2.2.9 was used for peak calling with default parameters. We detected a total of 20395 cells with 263346 non-overlapping peaks. Downstream analyses, such as quality control, filtering, normalization, integration, clustering and visualization was performed using the R package Seurat v5.0.170 and Signac v1.12.071. The appropriate QC matrices were calculated for each sample separately as per the standard Signac guidelines. The sample-specific putative doublets were predicted using AMULET v1.172 and filtered out from the subsequent analysis. We retained only high-quality cells (n=16546) based on number of fragments falling in peak regions (>3500), percentage of reads in peaks (>30%), and TSS enrichment (>2). The selected threshold was derived by manual inspection of the data distribution. To integrate the different samples an integration workflow was followed using high quality cells. Briefly, the four Prrx1 +/+ samples were merged using merge function provided by Seurat and then Term-Frequency Inverse-Document-Frequency (TFIDF) was calculated for the merged object. FindTopFeatures was used to select the top features with min.cutoff = 20. The Singular Value Decomposition (SVD) was performed on the selected top features. Harmony was used to perform the integration in SVD space. The first component was not considered for any downstream analysis as it is highly correlated with the sequencing depth in sATAC-seq data71.

Dimensionality reduction and clustering: A Shared Nearest Neighbor (SNN) Graph was built based on the harmony components using FindNeighbors function with top 25 components, and cells were clustered by a SNN modularity optimization-based clustering algorithm implemented in FindClusters function with resolution 0.15, this resulted in 6 clusters.

Imputation of gene expression and cluster annotation: To impute the gene expression matrix for sATAC-seq assay we employed the reference mapping approach described by73. First, the GeneActivity function was used to calculate gene activity score per cell in the sATAC-seq assay for the high variable genes used in scRNA-seq data. The gene activity score represents the count of fragments located within the 2 kb upstream region from the TSS of a gene. Next, the canonical correlation analysis (CCA) based anchors were identified between the scRNA-seq and sATAC-seq assay using FindTransferAnchors function. The scRNA-seq assay was used as a reference for transferring the gene expression matrix to sATAC-seq assay using TransferData function. The imputed gene expression values were visualized on UMAP using FeaturePlot function to show the specificity of the markers. Additionally, the genomic snapshots of open chromatin region around selected markers were visualized using CoveragePlot function.

Clustering of the cancer cell subset: To explore cancer cell heterogeneity at the open chromatin level, we subset cancer cells based on imputed tdTomato expression (also supported by other cancer cell markers). We re-integrated the data as explained in “Quality control, integration, dimensionality reduction and clustering” of sATAC-seq data analysis. Then, the neighbor searching was performed based on the top 20 harmony components and the FindClusters function was used with the resolution 1.0 to cluster the cancer cells into 16 distinct clusters. Next, we used similar gene signature-based approach as described in “Clustering of cancer cell subset” of scRNA-seq data analysis to annotate these clusters. The scores were calculated using imputed gene expression values for in sATAC-seq assay.

Motif enrichment analysis of the invasive cluster: The CORE collection of motifs for the vertebrates deposited in JASPIAR database was accessed using JASPAR2024 v0.99.6 R package74. The motif information was added to the Signac object using AddMotifs function and the motif activities was calculated per cell using chromVAR75 implemented in RunChromVAR function. Finally, the differential motif activity was calculated between invasive cluster (cluster15) and the rest of the cancer cell clusters using FindAllMarkers with mean.fxn = rowMeans and max.cells.per.ident = 500.

Prrx1 binding sites in genes of interest: To check whether Prrx1 has direct binding in the promoter regions of the candidate genes for dormancy, cyclins and cyclin inhibitors in the invasive cluster, first open accessible regions were identified in this cluster (cluster 15) using the AccessiblePeaks function. These accessible regions were then annotated to the nearest gene using the annotatePeak function from the ChIPseeker v1.38.0 package76. Promotor peaks associated with cyclins, cyclin inhibitors, and dormancy-related genes were analyzed for the presence of Prrx1 motif (JASPAR Database motif accession: MA0716.2), using the FIMO v5.5.777 motif search tool from the MEME suite. Coverage plots for these binding sites were generated using the CoveragePlot function.

In silico analysis of human cancer cell lines
The classification of 72 breast cancer cell lines along the epithelial (E) to mesenchymal (M) axis was taken from19. Briefly, the enrichment of epithelial (E) and mesenchymal (M) components were computed using epithelial and mesenchymal gene signatures27 in bulk gene expression data of breast cancer cell lines39. The EMT-TFs were excluded from the Mesenchymal signature to avoid biased correlations in the subsequent analyses. Additional scores such as, PING19, Proliferation (https://reactome.org/content/detail/R-HSA-69278), Dormancy35 and EMT (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp) signatures along with expression of cyclins, cyclin inhibitors, dormancy markers and PRRX1 were calculated in a similar fashion.

Cell Culture
BT-549 (ATCC) and BT-549 sh PRRX120 human tumor cell lines were cultured in DMEM (Sigma): F12 HAM (Sigma) media (1:1), supplemented with 10% heat-inactivated FBS (Sigma), 10 µg/ml insulin (Roche), 1% Penicillin-Streptomycin (Sigma), and 1% amphotericin (Sigma). MCDK-NBL2 cells (ATCC) were cultured in DMEM (Sigma) supplemented with 10% heat-inactivated FBS (Sigma), 1% Penicillin-Streptomycin (Sigma) and 1% amphotericin (Sigma). Cells were grown at 37°C and 5% CO₂, with media replaced every two to three days. Cells were passaged when they reached 80–90% confluency. For RNA extraction, 2x105 BT-549 cells were seeded in six-well plated in complete media and cultured for 48 hours and 105 MDCK-NBL2 cells were seeded in six-well plates in complete media and cultured for 72 hours.

Stable Prrx1 overexpression cell lines
Mouse Prrx1 was cloned into the pBABE-neo retroviral vector. Retroviral particles were generated in HEK-293T cells by transfection with Lipofectamine 2000 (Thermo Fisher Scientific) together with the packaging plasmids psPAX2 (RRID:Addgene_12260) and pMD2.G (RRID:Addgene_12259). Viral supernatants were collected 48 hours post-transfection and used to infect MDCK cells. Forty-eight hours after infection, MDCK cells expressing either pBABE-neo (empty vector) or pBABE-neo–Prrx1 were selected with G418 (400 µg/mL) for two weeks. Stable clonal cell lines were subsequently generated by limiting dilution and expanded for downstream analyses.

Transwell cell migration assay
The top chamber insert (Corning Costar Transwell) was covered with 50 μl of mouse collagen IV (Corning, 50 μg ml–1) and left to dry overnight. The resulting matrix was hydrated with 25 μl of water before seeding the cells. 5 × 104 cells from MDCK-NBL2 clones expressing different levels of Prrx1 were seeded and allowed to migrate. Cell nuclei at the bottom of the insert were imaged 8 hours after seeding and automatically counted using ImageJ.

Growth assay
A total of 5 × 10⁴ MDCK-NBL2 cells from clones expressing different levels of Prrx1 were seeded into 24-well plates in complete growth medium. Plates were placed in the Incucyte live-cell imaging system (Sartorius), and phase-contrast images were acquired every hour for 72 h. For each condition, three technical replicates were included in three independent experiments. Cell growth was quantified using the Incucyte’s built-in confluency analysis pipeline.

PRRX1 ChIP-Seq Data analysis
The raw data for PRRX1 ChIP-Seq in HS578, MDAMB157, and MDAMB436 cell lines were downloaded from Gene Expression Omnibus (GEO) database deposited under GSE202767 accession id. The raw data were filtered using Trimmomatic v0.39 tool with default parameters. The obtained high-quality reads were then aligned to the human reference genome (hg38) using bowtie v1.0.0 with default parameters. The aligned reads were sorted using samtools v1.11.0, and Picard MarkDuplicates v2.23.0 (https://broadinstitute.github.io/picard/) was used to identify and mark potential duplicated reads. To generate genomic tracks for PRRX1 ChIP-Seq signal for the selected genes deeptools bamCoverage v3.5.2 was used to convert the bam files into normalized coverage files (BigWig). Finally, the BigWig files were visualized in Integrated Genome Visualizer (IGV) v2.16.0 and genomic snapshots were taken for the selected genes.

Clinical analysis of gene signature in METABRIC data
The gene expression data for METABRIC cohort5,24,25 along with the clinical parameters was downloaded from cBioPortal data. The median of log normalized values was taken for the proliferation (Cell_Cycle_Mitotic_R-HSA-69278) and PINGs (Pro-Invasive genes) signatures. The patients were categorized based on quartiles, into high (75%) and low (25%) groups combining proliferation (R-HSA-69278) and invasion scores (PINGs). Finally, GraphPad Prism 9.4.1 was used to perform overall survival (OS) analysis and p value was calculated using the Gehan-Breslow-Wilcoxon test.

Statistics & Reproducibility
All experiments were repeated at least three times, and the nature and number of independent experimental replicates is indicated for each experiment in the figure legends. No data were excluded except in scRNA-seq, scATAC-seq and spatial transcriptomics analysis analyses. In those, bad quality cells were filtered out based on standard practice, as described in the Methods. No statistical methods were used to predetermine sample sizes. Sample sizes were empirically determined by similarity to those reported in previous publications from our laboratory using similar in vitro and in vivo experimental models. The experiments were not randomized. The investigators were not blinded to allocation during experiments and outcome assessment.
Statistical analyses were performed in R with RStudio, or in GraphPad Prism (GraphPad 9.4.1 software). Precise statistical tests are indicated in the figure legends. P values of >0.05 were considered not statistically significant (indicated as “ns”). Box plots are presented as means ± SEM. Boxplots show median, upper and lower quartile; whiskers show upper and lower extremes. For overall survival (OS) curves, p-values were obtained by using the Gehan-Breslow-Wilcoxon two-tailed test.

Data availability
The sequencing data generated in this study have been deposited in the Gene Expression Omnibus (GEO) repository. Source data are provided with this paper.
Spatial transcriptomics data that support the findings of this study have been deposited and are publicly accessible at the GEO repository:
(GSE291210, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE291210).
scRNA-seq data that support the findings of this study have been deposited and are publicly accessible at the GEO repository:
(GSE291211, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE291211).
snATAC-seq data that support the findings of this study have been deposited and are publicly accessible at the GEO repository:
(GSE291199, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE291199).

Code availability.
Detailed code pipeline has been deposited at a Zenodo repository:
https://doi.org/10.5281/zenodo.17726450

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Jimenez-Castano et al.\nThis is an interesting study that explores the role of Prrx1 transcription factor in regulation of the metastatic fate of\nbreast cancer cells. Using both human breast cancer and PyV mT GEMM, the authors propose that intermediate levels\nof Prrx1 are optimal to promote the metastatic colonization by balancing both invasive and proliferative potential of the\nemerging cancer cell. This is of potential interest; however the authors should clarify several issues.\n1. The authors use a K14-Cre reporter to excise Prrx1 function in mT tumor cells-As the tumors that are of luminal origin\nand should not express K14-The authors further claim that this reporter hits a stem cell population which should result in\nuniform labeling of all cancer cells. However, Fig 3d does not appear to show uniform red labeling. This raises the\nfurther question on what cell type is Prrx1 being deleted? The recent development of PyV mT IRES Cre GEMM where\nCre is coupled to PyV mT expression would perhaps address this issue.\n2. In a related issue dsRed is highly immunogenic. Can the authors preclude an immune impact on the observed\nmetastatic data.\n3. The premise of balanced proliferation and invasive programs ultimately determines metastatic fate. Have the\nauthors validated the contributions Prrx1 target genes are involved in this process.", "reviewer_2": "The study conducted by Jimenenz-Castano et al. contributes to the ongoing research by the Nieto group on the\nmechanisms underlying the impact of Prrx1 on metastatic dissemination. The study elucidates how Prrx1 functions in a\nhormetic manner, with intermediate expression levels associated with proliferating EMT phenotypes (ELM hybrids) and\nhigh levels reflecting enhanced dissemination of dormant DTCs. The authors employed a comprehensive approach that\nincorporated cell culture, mouse models, and human tissues to establish a correlation between varying Prrx1 levels and\nmetastatic outcomes. These findings are rigorous and well-controlled and generally support the primary tenets and\nconclusions of the study. A particularly intriguing finding that warrants further clinical validation pertains to patient survival\nand the association of varying Prrx1 expression with high and low scores for breast cancer cell proliferation, invasion, and\ndormancy signatures. While these and other findings are may contribute to the unification of several unresolved\nconundrums in the field, the study lacks novel mechanistic insights into the initiation, maintenance, and eventual\novercoming of hormetic Prrx1 expression. Furthermore, and based on an extensive scientific literature, several human cell\nlines identified as mesenchymal and dormant (Fig. 5g) exhibit full competence to complete the metastatic process in mice.\nConsequently, it presents a challenge to rectify hormetic Prrx1 function in these in silico analyses without additional\nvalidation experiments. Thus, it is advisable to undertake further attempts at validating Prrx1 levels in relation to their\nacquisition of proliferative, dormant, and invasive phenotypes in mice. Finally, monitoring alterations in Prrx1 splicing (e.g.,\nPrrx1a and Prrx1b) across the hormetic spectrum of Prrx1 function should be explored. Specific comments or concerns are\nprovided below under the section titled “Specific Comments.”\nSpecific Comments:\n1) Although the concept of hormetic signaling mediated by Prrx1 is intriguing, further development and empirical evidence\nare required to substantiate this claim. As previously mentioned, several of the cell lines depicted in Fig. 5g possess the full\ncapacity to complete the metastatic cascade in mice, including BT549 cells. Consequently, additional manipulation of Prrx1\nexpression in a dose-dependent manner is necessary to comprehensively assess the impact of low, intermediate, and high\nPrrx1 expression within the hybrid and mesenchymal groups. This is particularly pertinent for cells classified as hybrid, of\nwhich the vast majority appear to be devoid of Prrx1 expression. Therefore, generating and characterizing breast cancer\nvariants engineered to express graded levels of Prrx1 is warranted.\n2) The authors refer to Prrx1 as a single entity; however, Prrx1 is subject to alternative splicing, with Prrx1a governing\nmetastatic outgrowth (i.e., analogous to intermediate Prrx1 expression levels) and Prrx1b governing invasive epithelial-to-\nmesenchymal transition (EMT) phenotypes (i.e., analogous to high Prrx1 expression levels). Furthermore, simultaneous\ndepletion of both Prrx1a and 1b isoforms significantly reduces tumor initiation, growth, and suppresses the number and\nincidence of circulating tumor cells (CTCs) and metastases. These findings are remarkably similar to those presented in the\ncurrent study. Therefore, some attempt by the authors to monitor Prrx1 splicing in human and murine breast cancer samples\nis warranted. For instance, do Prrx1-high, mesenchymal breast cancer cells express robust levels of Prrx1b, while Prrx1-\nintermediate, hybrid breast cancer cells express robust levels of Prrx1a? Additionally, what is the impact of isoform-specific\ntargeting on breast metastasis, EMT, and dormancy programs? This is particularly important from a clinical standpoint, as\nisoform-specific detection may prove to be more effective in predicting disease trajectory and outcomes.\n3) A potentially intriguing finding that warrants further discussion pertains to the significant variations in stromal cell\nexpression of Prrx1 across the spectrum of its expression in carcinoma cells. Notably, Prrx1 has been identified as a CAF-\nspecific transcription factor, whose presence correlates with the development of chemoresistance and unfavorable disease\noutcomes. Consequently, it is pertinent to investigate whether differential expression of Prrx1 in CAFs compared to that in\ncarcinoma cells contributes to hormetic Prrx1 function and, consequently, disease outcomes. Some efforts to integrate these\nconcepts appear to be warranted.\n4) It is warranted to investigate molecular pathways that drive and maintain hormetic Prrx2 expression levels, such as\nTGFbeta and Wnt. Furthermore, it is crucial to determine whether targeted inactivation of these drivers impacts breast\ncancer dissemination and recurrence.\n5) Figure 5h should illustrate the extent of Prrx1 knockdown. Furthermore, similar augmentation and depletion of\nPrrx1 expression should be extended to additional hybrid and mesenchymal cell lines to demonstrate the plasticity\nassociated with varying levels of Prrx1 expression.", "reviewer_3": "This study addresses a valuable tumour biology question: how breast cancer metastatic potential is determined. Using a\nmulti-omics approach integrating patient cohorts, genetically engineered mouse models, and cutting-edge single-cell and\nspatial analysis technologies, the authors suggest that the transcription factor PRRX1 may play a crucial role as a key\nregulator of the invasion, proliferation, and dormancy programs. This study was generally well-designed, and the findings\nare novel, contributing to an improved understanding of tumour biology and potentially leading to therapeutic applications\nand the development of novel treatments. The authors highlight the importance of the hormetic transcriptional regulatory\nprogram in invasion, proliferation, and dormancy break. Moderate expression of PRRX1 confers the highest metastatic\npotential, enabling invasive cells to escape dormancy and proliferate. High expression of PRRX1 promotes invasion, but\nmay also enhance dormancy, potentially suppressing metastasis. Low or absent PRRX1 expression is associated with\nluminal subtype and reduced metastasis.\nThis study introduces the concept of hormetic regulation of metastasis. While hormetic regulation itself is not novel, its\nrole in invasion, proliferation, and dormancy breaks is intriguing. The hypothesis and mechanistic analysis suggest that\nseemingly contradictory biological activities interact in complex ways depending on the expression level of a single\nregulator, ultimately significantly affecting the behaviours of cells and cancer.\nMajor Concerns\nIn the MTA analysis, 94 cases of grade 3 IDC were selected. Selecting only grade 3 cases may drive bias in the\ninterpretation of results, at least from a clinical point of view. A larger cohort including grade 1 and 2 cases is needed.\nAdditionally, while stage 1, 2, and 3/4 cases were selected at roughly equal rates, the presence of stage 3/\n4\nmay influence\nthe analysis results. The size of the tumour burden is closely related to issues related to intra-tumoral heterogeneity and\nclonal evolution, which were carefully addressed in this study. Treatment guidelines may also have been different. It would\nbe desirable to describe the systemic therapy institutional concepts in 2010-2011. Additionally, the treatment details for the\nclinical cohort, particularly drug therapy, should be described.\nDifferential gene expression analysis between cancer cells in the invasive edge and bulk of the tumour is important, and the\nresults are interesting. Well-designed methods were used. This reviewer would like to know whether the same concepts and\ntechniques can be used to approach primary and metastatic tumour tissues. A similar question may be raised about spatial\ntranscriptome analysis, whether the same approach can be applied to other tissue microenvironments where dormancy\nbreaks occur.\nIs there any analysis on the two isoforms of PRRX1?\nWhat is the inter-observer concordance in pathologists' evaluations?\nRegarding the METABRCIS database analysis, what is the relationship of PRRX1 to major chromosomal abnormalities?\nWhat is the relationship with hormone-dependent growth mechanisms, including oestrogen receptor and FOXA1 axis?\nWhat is the relationship with pre- and post-operative adjuvant treatment? The relationship with recurrence and prognosis is\ncomplex, influenced by various factors, including tumour stage, biological characteristics, treatment sensitivity, and the host's\nimmune status. Therefore, jumping to conclusions about the relationship with clinical outcomes may have a significant impact\non hypothesis building and interpretation of results. Approaches such as examining a larger sample size by recurrence risk,\nstratifying treatments are necessary.\nRegarding hermetic role of Prrx1, is there a possibility that it may be modified by treatment?\nHow does the addition of cell cycle regulators, including CDK46 inhibitors, to the BT549 cell line (with or without hormone\ntherapy) affect these levels? This will likely be an important point in this study.\nIs it possible to directly confirm the regulatory role of PRRX1 through additional functional assays (e.g., knockdown or\nrescue experiments targeting dormancy-related genes such as GAS6/MME)?\nIntracellular metabolism, particularly cellular energy metabolism, appears to be important in the relationship between Prrx1\nexpression levels and EMT, invasion, proliferation, and dormancy. If authors have any knowledge on this point, it will\nstrengthen the hypothesis that the authors address.\nIn addition, how have the authors considered the interactions between tumour cells with different Prrx1 expression levels at\nthe invasion frontier? For example, what about energy exchange? What about interactions between tumour cells and stromal\ncells?\nFrom the perspective of Prrx1 expression inducers or inhibitors, what conditions would result in dormancy break? It will be\nbetter to explain, even as a hypothesis. The relationship between various treatments, such as endocrine therapy and anti-\nHER2 therapy, and dormancy break has been shown to be a long-term clinical phenomenon, in meta-analyses that have\ncompiled many prospective cases.\nThis may lead to the development of new treatments that take into account the regulation of PXRR1 expression.\nStatistical Robustness: As indicated, in some figures (e.g., metastasis incidence by PRRX1 level), the sample size is small.\nThe authors should provide confidence intervals etc.\nMinor Concerns\nRegarding the nonlinear correlation between PRRX1 expression level and metastasis incidence in breast cancer patients, the\npaper states that, while TNBC is known to include multiple subtypes, this data suggests that the high PRRX1 expression group\nmay specifically belong to a TNBC subtype with a low rate of metastasis. However, in Figure 1C, it is unclear how many TNBC\ncases actually existed among the basal cases with high PRRX1 expression (n=8) and which TNBC molecular subtype they\nbelonged to.\nIn the METABRIC database analysis, the negative correlation with growth appears weak. How definitive can this be?", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.", "reviewer_5": "In the manuscript “A hormetic transcriptional program coregulates invasion, proliferation and dormancy to define metastatic\npotential”, Jiménez-Castaño and colleagues investigate the metastatic mechanism in human breast cancer, specifically how\nexpression levels of the transcription factor PRRX1 regulate metastasis. To address the non-linear role of PRRX1 expression\nin correlation with human breast cancer incidence, the authors employed the MMTV-PyMT mouse model carrying floxed Prrx1\nalleles and performed spatial transcriptomics (MERFISH), single-cell RNA-seq, and single-cell ATAC-seq analyses. While the\noverall concept of investigating the role of PRRX1 in metastasis is interesting and of potential relevance to the field, the\nmanuscript remains largely descriptive with lack of in-depth bioinformatic analyses. Moreover, the limited functional and\northogonal validations further weakens the mechanistic insights, which limits its suitability for publication in Nature\nCommunications.\nAdditional points to be addressed:\n- Fig1 and Fig2, the use of IHC provides important insights into protein expression; however, the absence of RNA-level\nvalidation (e.g., RNAscope) limits the ability to determine whether the observed protein changes are transcriptionally\nregulated. Incorporating RNA-based analyses, such as single-cell or spatial transcriptomics, would significantly strengthen\nthe conclusions.\n- Extended Fig1, the authors should show each individual replicate for the cell boundary staining with DAPI, including\nzoomed-in views to demonstrate clear segregation of cells. In addition, quality statistics such as average detected genes\n(nGene) and counts (nCounts) for each replicate should be provided.\n- ATAC-seq quality controls, Key metrics are missing, including unique fragments per cell, fraction of reads in peaks (FRiP),\ntranscription start site (TSS) enrichment, fragment size distributions, and mitochondrial gene proportions. These are essential\nfor assessing dataset quality.\n- Fig 5a-f, when plotting the marker gene peaks, all populations should be visualized to show its specificity of the population of\ninterest.\n- Fig 5i, the ChIP-seq IGV tracks of different cell lines should be normalized and visualized on the same scale to allow direct\ncomparison.\n- Detailed Code should be provided.", "reviewer_6": "In this manuscript, et al. used a combination of human breast cancer data, genetic mouse models, and multi-omics data to\ncharacterise an invasive Prrx1+ state in primary breast cancer that has higher or lower metastatic potential based on the\nlevel of Prrx1 expression. They show that cells with an intermediate Prrx1 expression can maintain proliferation and an\nEMT-like invasive state, whereas high Prrx1 expression drives dormancy, showing that Prrx1 is an activator of cell cycle\ninhibitors and dormancy-associated genes.\nThe study provides very interesting and novel insights, supported by data of very good quality and strong validation.\nHowever, I have a few questions related mostly to the analyses of the data:\nFig3f - there are some stromal genes among the upregulated genes, such as Col1a1, Postn, Tnc - how sure are you that\nthose are cancer cells and not CAFs?\nFig 4 - Cluster 25 seems to be the one with the highest Prrx1 expression. Is this enriched in a specific genotype,\nfrom extended figure 3c it doesn’t see so, but can you show one umap coloured by genotypes, or 3 maps as in extended figure 3c\nbut coloured by cluster to make it clearer?\nFig 4 - Cluster 20 seems to have almost no Prrx1 expression, but still has high EMT and invasive score, how do you explain\nthis?\nFor the scRNA-seq analyses, the authors write that “The top 4000 HVGs were selected for the Prrx1 f/f and Prrx1 +/f\nsamples, whereas top 6000 HVGs were selected for Prrx1 +/+ samples”. Can they expand on why this was done?\nFig 5 - the authors imputed gene expression in scATAC-seq data using the scRNA-seq as a reference. Do they see a similar\nenrichment of the STINS score if they use the gene activity scores instead of imputation?\nMinor:\nAs someone who doesn't work with breast cancer, it was not clear why the authors used Krt14 areas at the tumour edge to\nidentify invasive areas, perhaps a very brief explanation could be added."}, "author_response": "Point-by-point response to reviewers\nReviewer 1\nWe thank the reviewer for the constructive review and for helping us to clarify some points. Here we\nprovide answers to the comments and suggestions.\n1. The authors use a K14-Cre reporter to excise Prrx1 function in mT tumor cells-As the tumors that are of\nluminal origin and should not express K14-The authors further claim that this reporter hits a stem cell\npopulation which should result in uniform labeling of all cancer cells. However, Fig 3d does not appear to\nshow uniform red labeling. This raises the further question on what cell type is Prrx1 being deleted? The\nrecent development of PyMT IRES Cre GEMM where Cre is coupled to PyV mT expression would perhaps\naddress this issue.\nWe appreciate the reviewer’s comment and apologize for not having described the mouse model in depth.\nThe reviewer refers to Fig. 3d and expected uniform labelling, given that the Cre recombinase is under the\ncontrol of the Krt14 promoter. However, the panel shown corresponds to Krt14 transcripts in advanced PyMT\ncarcinomas, restricted to basal-like cancer cells. This panel does not represent Tomato reporter or CRE\nexpression and we apologize if this was not clear in the previous version.\nThe key feature of the K14-Cre system is that Krt14 is expressed at the onset of embryonic mammary\ngland development, when mammary progenitors are uniformly K14+ (Wuidart et al., Nat Cell Biology 2018\n10.1038/s41556-018-0095-2). Cre activity at this stage activates the Tomato reporter and deletes Prrx1\nexon 2 in these progenitors, and the modification is subsequently inherited by all epithelial lineages,\nregardless of their later Krt14 expression status. For this reason, this lineage-wide recombination strategy\nhas allowed to study mammary gland development and stem cell contributions (Wuidart et al., Nat Cell\nBiology 2018 10.1038/s41556-018-0095-2) and we have successfully used it previously to label cancer\ncells in the PyMT model (Youssef et al., Nat Cancer 2024 10.1038/s43018024-00839-5).\nWe now include a new Extended Data Fig. 2 (see page 4 of main text) to better describe this Cre system\nin the PyMT breast cancer context, and showing the experimental design (a); (b) a schematic\nrepresentation of the temporal dynamics of K14-Cre activation and PyMT-driven transformation; (c)\nimmunofluorescence images of a normal adult mammary gland with all principal cells expressing\nTdTomato and specific epithelial markers (K8+K14- for luminal and K8-K14+ for myoepithelial cells).\nStromal cells are tdTomato-, confirming specific recombination in early mammary progenitors; (d) analysis\nof PyMT-driven carcinomas, showing that all PyMT+ cancer cells are Tomato+ and therefore bearing Prrx1\ndeletion since early development, and (e) quantification showing that nearly all K14+ cells (~99%) co-\nexpress PyMT, demonstrating that they derive from the transformed luminal-derived cancer population.\nWe believe that this discards the concern that K14+ cells in advanced carcinomas might represent\nuntransformed myoepithelial cells. It is also consistent with our recent scRNA-seq analysis using this\nmodel, and showing the reactivation of embryonic programs in PyMT tumors, derived from luminal cells\n(Youssef et al., Nat Cancer 2024 10.1038/s43018-024-00839-5). We hope this expanded explanation and\nthe new supplementary data clarify the model in which Prrx1 is efficiently and specifically deleted in\ncancer cells from the onset of tumor development.\n.\n2. In a related issue dsRed is highly immunogenic. Can the authors preclude an immune impact on the\nobserved metastatic data.\nWe thank the reviewer for highlighting the potential effects of fluorescent reporters on metastatic outcomes.\nWe agree that expression of exogenous proteins such as GFP or DsRed can increase immunogenicity,\npotentially reducing metastatic efficiency due to immune-mediated clearance (Grzelak et al., Cancer Cell 2022\n10.1016/j.ccell.2021.11.004). However, this effect is primarily observed in grafting or experimental metastasis\nassays, where cells expressing foreign antigens are introduced into mice with no prior exposure to those\nproteins.\nIn our system, tdTomato expression is activated early during embryonic development via the K14-Cre\nsystem and therefore being exposed to the immune system from the beginning. As previously mentioned,\nthis Cre system has been used to study normal mammary gland development without any reported\nimmunogenicity or developmental abnormalities (Wuidart et al., Nat Cell Biology 2018 10.1038/s41556-\n018-0095-2). Additionally, all experimental conditions in our study were conducted in the tdTomato\nbackground, ensuring that reporter expression is consistent across all groups and cannot account for\ndifferences in metastatic burden observed at varying Prrx1 levels.\nTo further address this concern, we analysed cytotoxic CD8+ T cell infiltration in metastatic sites where\nmetastases are consistently detectable (Prrx1 +/+ and Prrx1 +/f; see Figure below). We observed very\nlow overall infiltration, indicating minimal immune rejection of metastases. Importantly, there was no\nsignificant difference between the two genotypes, suggesting that the increased metastatic burden\nobserved in Prrx1 +/f tumors is not attributable to differential immune responses against tdTomato.\nImages and quantification of CD8 T cells (green) in lung metastasis (magenta) from Prrx1 +/+ and Prrx1 f/+ PyMT\nmouse breast tumors. Infiltrating CD8 cells (arrowheads) defined as green cells in contact with metastatic cancer\ncells (magenta). FC is represented as mean ± SEM. Prrx1 +/+ n= 8 metastases from 3 mice; Prrx1 +/f n= 7\nmetastases from 3 mice. Mann-Whitney test. ns, not significant. Scale bars: 50µm.\nIn summary, these results support that the metastatic phenotypes we observe are driven by Prrx1 dosage\nrather than by immunogenic effects of the reporter. Since this analysis serves as a methodological control\nrather than a result contributing to the main biological conclusions, we have included it here but not in the\nmain manuscript. We would be pleased to do it as an Extended Data Figure if the reviewer deems it\nnecessary.\n2\n3. The premise of balanced proliferation and invasive programs ultimately determines metastatic fate. Have\nthe authors validated the contributions Prrx1 target genes are involved in this process.\nThe Prrx1 target genes we describe in this manuscript are well-established regulators of proliferation\nand invasion. Cyclin D1/2 and the CDKN2 families are essential regulators of the G1/S transition.\nCCND1 is one of the most frequently amplified genes and CDKN2A/B is among the most frequently\ndeleted loci in breast cancer, see link below.\n(https://www.cbioportal.org/study/summary?id=brca_metabric).\nSimilarly, collagens and metalloproteases have long been validated as central mediators of invasive behavior.\nBecause the biological roles of these genes are well established, we believe that additional functional\nvalidation would not provide new insight into Prrx1-specific mechanisms. Regarding the dormancy-associated\ngenes, functional experiments were originally performed by the authors of the signature, please see the\nanswer to comment 11 by reviewers 3/4 for additional information.\nTo directly address the contribution of Prrx1 levels to invasion and cell-cycle control, we have now\nengineered a cell system expressing graded Prrx1 levels. These experiments demonstrate that differential\nPrrx1 expression modulates invasive capacity, proliferative arrest, and the expression of the proposed\nPrrx1 target genes (Main Figure 6d-g and discussed in page 12 main text), supporting the functional\nrelevance of the identified targets within the context of Prrx1-driven phenotypes.\nReviewer 2\nWe thank the reviewer for the appreciation of the different models used and for finding our analyses\nrigorous and well-controlled. Here we provide point-by-point answers to the comments and suggestions.\n1) Although the concept of hormetic signaling mediated by Prrx1 is intriguing, further development and\nempirical evidence are required to substantiate this claim. As previously mentioned, several of the cell\nlines depicted in Fig. 5g possess the full capacity to complete the metastatic cascade in mice, including\nBT549 cells. Consequently, additional manipulation of Prrx1 expression in a dose-dependent manner is\nnecessary to comprehensively assess the impact of low, intermediate, and high Prrx1 expression within\nthe hybrid and mesenchymal groups. This is particularly pertinent for cells classified as hybrid, of\nwhich the vast majority appear to be devoid of Prrx1 expression. Therefore, generating and characterizing breast\ncancer variants engineered to express graded levels of Prrx1 is warranted.\nWe thank the reviewer for raising this important point regarding PRRX1 hormetic signaling and the\ninterpretation of human breast cancer cell line data.\nThe reviewer specifically refers to Fig. 5g (now Fig. 6a) and mentions BT549 cells, noting their reported\nmetastatic capacity. We understand that laboratory conditions, culture history, and potential genetic\ndrift can influence outcomes in experimental metastasis assays. However, we have previously shown\nthat parental BT549 cells are unable to form metastases in vivo, and metastatic behaviour is only\nobserved upon experimental reduction of Prrx1 (Ocaña et al., Cancer Cell 2012\n10.1016/j.ccr.2012.10.012, see Figure below).\nRepresentative GFP-positive foci in lungsand H&E sections of the lungs. Scale bars: 1 mm (upper panels)\nand 200 μm (lower panels).\nThis is further corroborated by the METMAP dataset (Jin et al., Nature 2020 10.1038/s41586-020-29692),\nwhich systematically characterizes hundreds of human cancer cell lines. BT549 cells display 0\npenetrance, indicating that barcodes from these cells could not be recovered in metastasis from any\nanimal, reflecting their very low intrinsic metastatic potential (see Figure below). Thus, we believe that\nBT549 remains an appropriate model to study the role of Prrx1 in regulating metastatic capacity.\na, Scheme of METMAP system to systematically score metastatic potential of human cancer cell lines (upper\npanel). b, Metastatic score, potential and penetrance, of the BT549 human breast cancer cell line. Extracted\nfrom the METMAP explorer and database: https://depmap.org/metmap/.\nthis is what we have done in this study by genetically engineering graded Prrx1 levels specifically in cancer\ncells in vivo in a spontaneous, immunocompetent mouse model of breast cancer, preserving the natural\nintratumoral heterogeneity, physiological Prrx1 dosage, and metastatic progression. This experimental design\nhas enabled us to perform a direct and physiological analysis of the impact of negative, intermediate, and\nhigh Prrx1 levels on metastatic capacity. Nevertheless, we agree that engineering cells to express graded\nlevels of Prrx1 is a pertinent addition to our data and we have developed a cell culture system that allows\ncontrolled assessment of low/negative, intermediate, and high Prrx1 expression states starting from an\nepithelial state using MDCK cells. We have decided to use this cell line because it is the most used cell line in\nEMT and plasticity studies in general and also to avoid peculiarities or aberrant responses typically associated\nwith individual cancer cell lines.\nThese new data are shown in revised Fig. 6d-g (see page 12 in main text). We show that:\n• Prrx1-negative cells proliferate and lack invasive properties.\n• Intermediate Prrx1 levels promote invasion without impacting proliferation, thereby maximizing\nmetastatic potential.\n• High Prrx1 levels induce a similar invasive behavior but strongly suppress proliferation and induce\ndormancy.\nAltogether, these experiments are fully compatible with the observed behavior associated with Prrx1 levels in\nour in vivo model.\nWe also acknowledge the reviewer’s comment on the visualization of Prrx1 levels across EMT groups. The\noriginal heatmap normalization and color scale may have obscured intermediate Prrx1 levels. We have now\nincluded a new plot for this data in Extended Data Figure 11a showing that as a group, the hybrid cells have\nPRRX1 transcript levels in between epithelial and mesenchymal populations.\nImportantly, the absence of Prrx1 in some hybrid cell lines does not contradict our conclusions. The\nhormetic effect of intermediate Prrx1 expression on metastatic potential is only revealed in in vivo\nmodels, where cells pass the full through the whole metastatic cascade from invasion to colonization.\nIn contrast, in experimental metastasis assays (for example tail-vein injections) cells bypass the early\nsteps, invasion and intravasation. In such assays, cell lines with low or no Prrx1 expression may still\ncolonize secondary sites efficiently.\n2) The authors refer to Prrx1 as a single entity; however, Prrx1 is subject to alternative splicing, with\nPrrx1a governing metastatic outgrowth (i.e., analogous to intermediate Prrx1 expression levels) and\nPrrx1b governing invasive epithelial-to-mesenchymal transition (EMT) phenotypes (i.e., analogous to\nhigh Prrx1 expression levels). Furthermore, simultaneous depletion of both Prrx1a and 1b isoforms\nsignificantly reduces tumor initiation, growth, and suppresses the number and incidence of circulating\ntumor cells (CTCs) and metastases. These findings are remarkably similar to those presented in the\ncurrent study. Therefore, some attempt by the authors to monitor Prrx1 splicing in human and murine\nbreast cancer samples is warranted. For instance, do Prrx1-high, mesenchymal breast cancer cells\nexpress robust levels of Prrx1b, while Prrx1-intermediate, hybrid breast cancer cells express robust\nlevels of Prrx1a? Additionally, what is the impact of isoform-specific targeting on breast metastasis, EMT, and\ndormancy programs? This is particularly important from a clinical standpoint, as\nisoform-specific detection may prove to be more effective in predicting disease trajectory and outcomes.\nWe thank the reviewer for raising the point of Prrx1 isoforms and their potential functional relevance. Indeed,\nPrrx1 undergoes alternative splicing of exons 3–5, giving rise to two isoforms, which we call Prrx1 Long (L)\nand Prrx1 Short (S) to avoid confusion with the two gene duplicates of Prrx1 (prrx1a and prrx1b) in teleosts.\nThe L isoform (ortholog to prrx1a) arises from alternative splicing between exons 3 and 5, producing a protein\nthat includes the AD3 and OAR domains, whereas the S isoform (ortholog to prrx1b) results from normal\nsplicing between exons 3 and 4, generating a shorter protein lacking these domains due to a stop codon in\nexon 4. See Figure below for a scheme of the Prrx1 loci and isoform structure.\na, Schematic representation of the Prrx1 gene, which has 5 exons. Alternative splicing of exon 4 results in the\ngeneration of two\ndifferent protein\nproducts. Noted that\nour genetic KO\ndeletes exon 2, which\nis common to both\nisoforms and contain\nthe critical homeobox\ndomain b, Schematic representation of the two proteins produced by the Prrx1 gene. The main difference\nbetween both proteins is that PRRX1L includes the AD3 and OAR domains, and PRRX1S does not.\nAware of this regulatory layer, our Prrx1 conditional knockout (cKO) design targets exon 2, which is shared\nby both isoforms and encodes the critical homeobox domain. We have previously demonstrated that this\nstrategy effectively abolishes both isoforms at the protein level in embryonic fibroblasts obtained from mouse\nembryos carrying this genetic design. Furthermore, heterozygous deletion reduces the levels of both isoforms\nin a comparable manner (Youssef et al., Nat Cancer 2024 10.1038/s43018024-00839-5) (see Figure below).\nPrrx1 expression in fibroblasts obtained from mouse embryos of different genotypes. Note that TAT-CRE\ntreated Prrx1 flox/flox embryos (well #3) has the same profile as Prrx1 homozygous mutant embryos (well #5),\nshowing one band that corresponds to the product of gene duplicate Prrx2. Importantly, both Prrx1 isoforms are\nabsent at protein level. Notice how Prrx1 heterozygous (flox/+) shows a similar reduction in both isoforms, with\nintermediate levels between the WT and the homozygous, as expected. The WB corresponds to MEFs obtained\nfrom one of three independent validation experiments. Extracted from : 10.1038/s43018-024-00839-5.\nbalance between isoforms, we designed isoform-specific primers and reanalysed FACS-sorted cancer cells\nfrom our experimental groups. This confirmed that both isoforms are reduced in a similar manner, and\nrelative to the number of available alleles ruling out differential isoform effects as the origin of the\nmetastatic phenotypes. Please see revised Figure 2a and page 5 of main text.\nTo specifically address the reviewer´s question, we measured isoform-specific expression by qPCR in eight\ncell lines spanning epithelial, hybrid, and mesenchymal groups based on our in-silico analysis. We did not find\nevidence of preferential isoform expression, and found that levels of both isoforms similarly increase with\nmesenchymalization (see Figure below, panel a).\nAdditionally, in a recent study, we treated MCDK NBL2 epithelial cells with TGFβ to induce EMT. Both isoforms\nincreased concomitantly with acquisition of invasive phenotypes (Youssef et al., Nat Cancer 2024\n10.1038/s43018-024-00839-5) (see Figure below, panel b). These data indicate that, unlike in pancreas, Prrx1\nisoforms do not have differential activation in our breast cancer or kidney models.\nIn summary, while isoform-specific regulation is an important consideration, our mouse and human breast\ncancer data indicate that both Prrx1 isoforms are coordinately expressed in breast cancer, and the critical\ndeterminant of metastatic behavior is the level of Prrx1, rather than isoform composition.\na, Relative transcript levels (RT-qPCR) for PRRX1 long (PPRX1-L) and short (PPRX1-S) in human breast cancer cell line along\nthe ENT spectrum. qPCR fold-change expression is normalized to NDA 436 levels and represented as mean ± SEN, n = 3.\nb, Relative transcript levels (RT-qPCR) for PRRX1 long (PPRX1-L), short (PPRX1-S) or both isoforms (PRRX1) in\nNDCK-NBL2 cells, either 1 or 4 days after TGFβ treatment. qPCR fold-change expression is represented as mean ± SEN,\nn = 3 independent experiments/time point and condition. Extracted from 10.1038/s43018-024-00839-5.\n3) A potentially intriguing finding that warrants further discussion pertains to the significant variations\nin stromal cell expression of Prrx1 across the spectrum of its expression in carcinoma cells. Notably,\nPrrx1 has been identified as a CAF-specific transcription factor, whose presence correlates with the\ndevelopment of chemoresistance and unfavorable disease\noutcomes. Consequently, it is pertinent to\ninvestigate whether differential expression of Prrx1 in CAFs compared to that in carcinoma cells\n7\ncontributes to hormetic Prrx1 function and, consequently, disease outcomes. Some efforts to integrate these\nconcepts appear to be warranted.\nWe thank the reviewer for considering the role of Prrx1 in other cell types. However, it is important to\nemphasize that our system is cancer-cell-specific (please see the new extended description of the mouse\nmodel in Extended Data Figure 2 and page 5 of main text). Our system does not target other cell types,\nincluding cancer-associated fibroblasts (CAFs), and the phenotypes we report are exclusively to Prrx1 deletion\nin carcinoma cells.\nThe reviewer rightly alludes to Prrx1 expression in fibroblasts, and indeed, it is also found in multiple\nother cell types, including pericytes (Pei et al., Nature Genetics 2025 10.1038/s41588-025-02341-9),\nmelanoma cells (Karras et al., Nature 2022 10.1038/s41586-022-05242-7) (Ferreres et al., Molecular\nOncology 2024 10.1002/1878-0261.13688), and pancreatic cancer cells (Takano et al, Genes &\nDevelopment 2016 10.1101/gad.263327.115). We acknowledge that, given our observations on Prrx1\nfunction in cancer cells, exploring Prrx1 levels in other cell types could be of great interest. However,\nsuch studies are beyond the scope of the current manuscript.\n4) It is warranted to investigate molecular pathways that drive and maintain hormetic Prrx2 expression levels,\nsuch as TGFbeta and Wnt. Furthermore, it is crucial to determine whether targeted inactivation of these\ndrivers impacts breast cancer dissemination and recurrence.\nWe appreciate the reviewer’s interest in understanding the upstream pathways regulating Prrx1. TGF-13\nand Wnt signaling have been extensively linked to EMT, and previous studies have shown that TGF-13\ninduces PRRX1 expression and that Wnt contributes to PRRX1-mediated EMT induction in the MDCK\nmodel (Youssef et al., Nat Cancer 2024 10.1038/s43018-024-00839-5). However, experimentally\ninactivating these pathways would affect hundreds of downstream processes and broad gene-regulatory\nnetworks, and in this study we focus on Prrx1 impact on the regulation of metastatic colonization through\nthe characterization of its downstream targets.\n5) Figure 5h should illustrate the extent of Prrx1 knockdown. Furthermore, similar augmentation and\ndepletion of Prrx1 expression should be extended to additional hybrid and mesenchymal cell lines to\ndemonstrate the plasticity associated with varying levels of Prrx1 expression.\nThanks for the suggestion. We have now included the quantification of Prrx1 knockdown in Figure 6b. In\naddition, to further support the plasticity and regulatory effects associated with varying Prrx1 levels, we\nanalysed publicly available data in HS578T cells, a mesenchymal breast cancer cell line with high endogenous\nPrrx1 expression (Figure 6a), from the same study that provided the Prrx1 ChIP-seq dataset (Jovanović et al.,\nCell Reports 2023 10.1016/j.celrep.2023.113564). The dataset includes bulk RNA-seq following Prrx1\ndepletion, and we found that among the statistically significant differentially expressed genes, several of them\noverlapped with those examined in our manuscript. Notably, CCND2, the cyclin D family member, was\nsignificantly upregulated upon Prrx1 knockdown, whereas CCND1 and CCND3 changed to a lesser extent.\nFurthermore, MME, one of the dormancy-associated markers downregulated after Prrx1 loss in vivo and in\nBT549 cells, was also significantly reduced (new Extended Data Figure 11b and page 10). This, together with\nthe new graded Prrx1 overexpression experiments in MDCK cells (Fig. 6d-g), provide further support to the\nmultiregulatory role of Prrx1 and the model presented in this study.\nWe thank the two reviewers for the very positive comments on the question addressed in this work, the\nexperimental approaches used, the novelty and potential therapeutic application. We also acknowledge the\nthorough review of the clinical data we include, and we very much appreciate the multiple questions and\nsuggestions with respect to those data. We have carefully read and addressed all the questions and we\nrespectfully would like to argue that many of them lie in the field of clinical research, which is far from the\nfocus of this work.\nMajor Concerns\n1) In the MTA analysis, 94 cases of grade 3 IDC were selected. Selecting only grade 3 cases may drive bias\nin the interpretation of results, at least from a clinical point of view. A larger cohort including grade 1 and 2\ncases is needed.\nPRRX1 expression has previously been associated with high-grade breast tumors (Ocaña et al., Cancer\nCell 2012 10.1016/j.ccr.2012.10.012). More recently, we demonstrated that PRRX1 expression is also\nlinked to triple-negative breast cancer (TNBC) (Youssef et al., Nat Cancer 2024 10.1038/s43018-024-\n00839-5), a subtype that is typically diagnosed as high grade and is strongly associated with tumor\nprogression and poor prognosis. Based on these findings, we aimed to further investigate the role of\nPRRX1 across molecular breast cancer subtypes that share similar histological characteristics and tumor\ngrade.\nIn addition, our previous studies have shown that PRRX1 expression is associated with metastatic behavior.\nThis association would be difficult to robustly assess in low-grade tumors, where aggressive subtypes are\nrarely represented. Therefore, restricting the analysis to grade 3 tumors allowed us to focus on a clinically\nrelevant population in which the biological effects of PRRX1 are more likely to be observed.\nFinally, the selected cohort includes extensive and well-annotated clinical information, which is essential\nfor evaluating the prognostic value of the analyzed marker. Taken together, these considerations support\nthe rationale for focusing on grade 3 IDC cases while minimizing confounding factors related to tumor\ngrade.\nAdditionally, while stage 1, 2, and 3/4 cases were selected at roughly equal rates, the presence of stage 3/4\nmay influence the analysis results.\nWe respectfully disagree with the assertion that the inclusion of stage III/IV tumors introduces bias into the\nanalysis. On the contrary, given that the cohort is intentionally restricted to grade 3 tumors, the inclusion of\ncases spanning stages I to IV was a deliberate strategy to reduce potential bias and to ensure adequate\nclinical representation across the full spectrum of disease progression.\nExcluding advanced-stage (III/IV) tumors would result in an incomplete and clinically unbalanced series,\nparticularly for high-grade breast cancers, which frequently present more advanced stages. Therefore, the\npresence of stage III/IV cases does not distort the analysis; we believe that it rather strengthens the study\nby allowing the assessment of PRRX1 expression within a biologically and clinically coherent high-grade\npopulation across different disease stages and subtypes.\nlow-grade breast cancers are managed in the neoadjuvant setting, and diagnostic material is often restricted\nto fine-needle aspiration (PAAF) or small biopsies, which are insufficient for the analyses performed here.\nIn addition, PRRX1 expression has been consistently associated with aggressive, high-grade tumors, supporting\nour focus on grade 3 disease.\nTaken together, the study design was carefully chosen to balance biological relevance and clinical\nheterogeneity. While we acknowledge that larger and more heterogeneous cohorts will be valuable in future\nstudies, we believe that the current cohort composition is appropriate and does not introduce any significant\nbias.\n2) The size of the tumour burden is closely related to issues related to intra-tumoral heterogeneity and\nclonal evolution, which were carefully addressed in this study.\nWe acknowledge that tumor burden can be associated with intra-tumoral heterogeneity and clonal evolution.\nHowever, we believe that this concern does not substantially affect the present study for several reasons\nrelated to breast cancer biology and to the pathological methodology applied.\nFirst, although intra-tumoral heterogeneity is a recognized phenomenon in cancer, breast cancer is not\ntypically characterized by marked genetic heterogeneity within primary tumors, particularly when\ncompared with other solid malignancies. While epigenetic heterogeneity may exist, immunohistochemical\nanalyses are routinely and reliably used in breast cancer diagnostics and research, and their\ninterpretability is well established in this context.\nSecond, potential intra-tumoral heterogeneity was explicitly addressed at the technical level. All tumors\nincluded in the tissue microarrays (TMAs) were sampled using two independent cores obtained from different\ntumor areas, following standard pathological practice. This approach is specifically designed to account for\nspatial heterogeneity and is routinely employed in breast cancer TMA-based studies. Heterogeneity related to\ntumor size or burden was not overlooked, but rather systematically considered through appropriate sampling\nand established immunohistochemical methodology. We believe that this strategy adequately mitigates the\npotential impact of intra-tumoral heterogeneity on the analysis and supports the robustness of the conclusions\ndrawn.\nTreatment guidelines may also have been different.\nWe agree that treatment guidelines differ according to breast cancer molecular subtype. Luminal tumors\nare typically treated with chemotherapy and/or endocrine therapy, HER2-positive tumors with\nchemotherapy combined with anti-HER2 agents, and triple-negative or basal-like tumors primarily with\nchemotherapy. These differences are well established and were taken into account when interpreting\nthe results.\nImportantly, all tumors included in this study were retrospectively collected in the adjuvant setting according\nto the corresponding guidelines at the time. The fact that PRRX1 expression is associated with invasive and\naggressive features across molecular subtypes, despite these distinct therapeutic regimens, supports the\nnotion that PRRX1 marks tumor biology related to invasion rather than reflecting treatment-specific effects.\nFurthermore, the inclusion of grade 1–2 tumors would not be comparable\n10\nwithin this framework. In current clinical practice, many low-grade breast cancers are managed in the\nneoadjuvant setting, with systemic treatment administered prior to surgery. As a result, post-treatment\nsurgical specimens from grade 1–2 tumors are often unavailable or biologically altered by therapy, precluding\nmeaningful comparison with adjuvant-treated grade 3 tumors. Therefore, the homogeneity of the cohort in\nterms of treatment timing (adjuvant setting) and tumor grade strengthens the internal consistency of the\nstudy and allows assessment of PRRX1 correlation with metastasis incidence independently of treatment\nmodality.\nIt would be desirable to describe the systemic therapy institutional concepts in 2010-2011. Additionally, the\ntreatment details for the clinical cohort, particularly drug therapy, should be described.\nWe agree that providing information on institutional systemic therapy concepts may be useful for context,\nand we have therefore clarified that patients were treated according to contemporaneous national and\ninternational clinical guidelines. However, it is important to note that the cohort includes patients\ndiagnosed over a broader time period, including cases from 2003–2004, and not exclusively from 2010–\n2011. While systemic treatment strategies have evolved over time, the backbone of chemotherapy\nregimens during this period, primarily anthracycline and taxane-based combinations, remained largely\nconsistent. As such, these variations are unlikely to substantially influence the biological associations\ninvestigated in this study.\nImportantly, and as discussed later in our response, mechanistic cancer-biology studies published\nrecently with similar experimental design (Rossi et al., Nature 2022 10.1038/s41586-022-04758-2)\ntypically do not include detailed systemic-treatment information when such data are not available or\nare not essential for the experimental and mechanistic focus of the study. Our manuscript follows this\nestablished practice.\nTo address the reviewer’s concern as rigorously as possible, we have now performed two analyses using the\nMETABRIC dataset. These analyses show that neither PRRX1 expression (see Figure below) nor our patient\nstratification based on the proposed signatures (Extended Data Figure 13) differ significantly across treatment\n. groups For these reasons, while we acknowledge the reviewer’s suggestion, we believe that a detailed\ndescription of year-specific institutional treatment algorithms would not alter the interpretation of the results.\nPRRX1 expression levels across different treatment groups in METABRIC.\n11\n3) Differential gene expression analysis between cancer cells in the invasive edge and bulk of the tumour is\nimportant, and the results are interesting. Well-designed methods were used. This reviewer would like to\nknow whether the same concepts and techniques can be used to approach primary and metastatic tumour\ntissues. A similar question may be raised about spatial transcriptome analysis, whether the same approach\ncan be applied to other tissue microenvironments where dormancy breaks occur.\nWe fully agree with the reviewer that the conceptual framework and methodological strategy applied here,\nleveraging spatially resolved transcriptomics to dissect region-specific cancer cell states, can be extended to\nother pathological contexts, including both primary vs. metastatic lesions and microenvironments where\ndormancy reactivation occurs. Spatially guided identification of invasive, dormant, or reactivating niches\nrepresents an exceptionally powerful approach to understanding metastatic progression. However, the aim\nof the present study was to specifically investigate how preexisting heterogeneity within the primary tumor\nrelates to metastatic potential. For this reason, all spatial and transcriptional analyses were focused on the\nprimary tumor, complemented by targeted validation of cancer-cell behavior at metastatic sites (as in the\ncase of metastatic mitotic ratio). Applying the same high-resolution spatial strategies to metastatic tissues\nand dormancy-break microenvironments is indeed of great interest and represents a highly promising\ndirection for future work.\nIs there any analysis on the two isoforms of PRRX1?\nWe have incorporated new analysis of PRRX1 isoforms in the manuscript (Figure 2a) and provide\nadditional information in this letter. Please check our answer to comment 2 of reviewer 2.\n4) What is the inter-observer concordance in pathologists' evaluations?\nInter-observer concordance between pathologists was high, exceeding 90%, in line with standard practice\nfor immunohistochemical evaluation in breast pathology. All cases were independently assessed by two\nexperienced pathologists. In instances where substantial discordance was observed, the immunostaining\nwas repeated and evaluated on full tissue sections rather than on TMA cores, allowing for comprehensive\nreassessment and consensus scoring. This approach ensures robustness and minimizes the impact of\nsubjective variability. In addition, for the classification based on PRRX1 (immunofluorescence) expression\npattern, we incorporated an external, independent validation using a digital pathology–based approach in\na subset of the analysed samples. This digital assessment showed 100% concordance with the\npathologists’ classification, further supporting the robustness of our scoring strategy. These procedures\nare routinely applied in our laboratory and are consistent with established quality control practices in\ndiagnostic and research pathology. This information has now been clarified in the Methods section of the\nrevised manuscript (page 27 and 29 of the Ms).\n6) Regarding the METABRCIS database analysis, what is the relationship of PRRX1 to major chromosomal abnormalities?\nWe recognize the importance of understanding the genomic landscape of PRRX1 in human breast cancer, as\nit may contribute to heterogeneity in expression levels. Motivated by the reviewer’s question, we performed\na deeper analysis of the METABRIC dataset that has been incorporated in the revised version in the manuscript\nas new Extended Data Figure 1 and discussed in page 3 and 4 of main text. In summary, we found that\nPRRX1 is one of the most frequently amplified genes in the METABRIC dataset, with amplification present in\n~20% of patients. Many of these patients are characterized by 1q gain and 16q loss, where the PRRX1 locus\nis located (1q). However, this genomic alteration encompasses\n12\nhundreds of genes, making it difficult to isolate the specific contribution of PRRX1 amplification to clinical\nor molecular features. In addition, this amplification does not seem to be translated to significantly higher\nlevels of PRRX1 transcripts.\n6) What is the relationship with hormone-dependent growth mechanisms, including oestrogen receptor and\nFOXA1 axis?\nPRRX1 is not a canonical component of hormone-dependent growth pathways and is not directly regulated\nby the oestrogen receptor (ER)–FOXA1 transcriptional axis. Instead, PRRX1 has been consistently linked\nto metastatic competence, largely independent of hormone-driven proliferation. In line with this, PRRX1\nexpression in our cohort was predominantly associated with ER-negative and basal tumors, which are\ncharacterized by limited reliance on oestrogen-dependent signaling and frequently exhibit EMT-related\nfeatures. This observation is consistent with previous studies demonstrating an inverse relationship\nbetween PRRX1 expression and luminal differentiation programs cells (Youssef et al., Nat Cancer 2024\n10.1038/s43018-024-00839-5).\nTo further address the reviewer’s point, we performed additional analyses of the METABRIC dataset,\nfocusing specifically on PRRX1, ESR1, and FOXA1 (see figure below). As expected, ESR1 and FOXA1 exhibit\na strong positive correlation, reflecting their known co-regulatory relationship. In contrast, PRRX1 shows\na modest but consistent negative correlation with both ESR1 and FOXA1, reinforcing the notion that high\nPRRX1 expression aligns with more dedifferentiated, non-luminal tumor states. Because ESR1 and FOXA1\ndisplay clear bimodal distributions across patient samples, direct linear correlation metrics have limited\ninterpretability. However, when correlations are examined within the ESR1-high / FOXA1-high subgroup,\nthe negative association with PRRX1 seems to be more pronounced, particularly for ESR1.\nTaken together, our data support the interpretation that PRRX1 marks a hormone-independent, invasive\ntumor state rather than modulating oestrogen-dependent growth mechanisms.\nCorrelation analysis between the expression of ESR1, FOXA1 and PRRX1. Selected from METABRIC breast\ncancer human data (n=1980). Red line represents regression line. Two-tale Pearson correlation test. **** p\n< 0.0001.\n13\n7) What is the relationship with pre- and post-operative adjuvant treatment?\nAll cases included in this study belong to an adjuvant-treated cohort. Tumor samples were obtained exclusively\nfrom surgical specimens collected after completion of standard-of-care adjuvant therapy. The study was not\ndesigned to compare PRRX1 expression in matched pre-treatment (diagnostic biopsy) and post-treatment\n(surgical) samples.\nWhile such longitudinal analyses would indeed be of interest, they are not feasible in the present cohort.\nDiagnostic biopsies are not available for all patients, and in some cases, patients have either deceased or are\nno longer followed at our institution, precluding systematic paired analyses. As a result, any pre-/post-\ntreatment comparison would be incomplete and potentially biased. We agree that future studies specifically\ndesigned to analyze matched pre- and post-treatment samples will be necessary to address dynamic changes\nin PRRX1 expression. However, such analyses would require the establishment of a new, prospectively\ncollected cohort.\n8) The relationship with recurrence and prognosis is\ncomplex, influenced by various factors, including tumour\nstage, biological characteristics, treatment sensitivity, and the host's immune status.\nWe agree that recurrence and prognosis in breast cancer are multifactorial and influenced by numerous\nvariables, including tumor stage, biological subtype, treatment response, and host-related factors. However,\nour analyses focused on established and clinically relevant tumor-related variables for which reliable data\nwere available. While additional patient-related factors such as body mass index, smoking status, menopausal\nstatus, or immune parameters may contribute to prognosis, incorporating all potential confounders would\nrequire a specifically designed, large-scale, multicenter prospectively study, which goes far beyond the scope\nof this work.\nImportantly, the consistent association of PRRX1 expression with invasive features across the analyzed tumor\nseries supports its relevance as a biological marker, independent of the broader array of factors that influence\npatient outcome. This limitation and the need for future prospective studies have now been acknowledged in\nthe revised version (pages 13-14).\nTherefore, jumping to conclusions about the relationship with clinical outcomes may have a significant impact\non hypothesis building and interpretation of results.\nWe fully agree that overinterpretation of clinical outcome associations should be avoided. Importantly,\nthis was not the intention of the present study. The manuscript does not aim to establish PRRX1 as a\nclinical decision-making or predictive biomarker, nor to draw causal conclusions regarding patient\noutcomes. Rather, the objective of the section of this work referring to human data and samples is to\nassess compatibility with the new biological data and concepts put forward in the mouse model, a factor\nthat collectively regulates cancer cells features that are key for metastatic potential, and all of this carried\nout in a well-characterized tumor series from patients. Any associations with recurrence or prognosis are\npresented in a descriptive and hypothesis-generating context, consistent with the study design and sample\nsize. As mentioned above, we acknowledge that definitive conclusions regarding clinical outcomes would\nrequire large, prospective, multicenter studies specifically designed to address treatment response and\nlong-term survival. This limitation has been clarified in the revised Discussion to prevent misinterpretation\nof the results (pages 13-14).\nnecessary.\nWe agree that analyses based on larger sample sizes, stratification by recurrence risk, and treatment-\nspecific subgroups would be valuable and would further strengthen the clinical interpretation of PRRX1\nexpression. However, such approaches are beyond the scope of the present study, which was designed\nto validate a biological marker in a well-annotated tumor series, not to perform treatment-stratified or\nlongitudinal outcome modeling. Implementing the analyses suggested by the reviewer would require the\nestablishment of independent validation cohorts, ideally including paired pre- and post-treatment samples\nand both neoadjuvant and adjuvant settings. This process involves not only access to biobank material,\nbut also rigorous selection of matched specimens, detailed clinical annotation, and extensive analytical\nwork, all of which require substantial time and resources.\nWhile we appreciate that the dataset may appear limited in comparison with dedicated clinical studies,\nthe scope and sample size are consistent with recently published high-impact studies. For example, in\n(Rossi et al., Nature 2022 10.1038/s41586-022-04758-2), the authors use a cohort of 129 breast cancer\npatients, comparable to our 94, and classify them into high, intermediate, and low PGDH expression\ngroups by IHC, without stratification by treatment and even without distinguishing the cell type of origin\nof the signal. That study similarly used this initial observation as motivation to investigate the mechanistic\nrole of PGDH in metastasis in vivo.\nIn addition, we have previously published correlations between PRRX1 and disease-free survival in both\nbreast and lung cancer using bigger multiple independent cohorts (Ocaña et al., Cancer Cell 2012\n10.1016/j.ccr.2012.10.012) (see figure in next page). However, those analyses were based on bulk RNA\nexpression, which cannot differentiate cancer-cell–specific expression and only consider two groups based\non Prrx1 expression. For this reason, we present here a refined analysis based on protein expression and\ncancer-cell–specific scoring, which provides a more accurate framework for interpreting the biological\nfindings.\nNevertheless, to address the reviewer’s concern as rigorously as possible, we have now incorporated two\nanalyses using the METABRIC dataset. These analyses show that neither PRRX1 expression (see Reviewer 2-\nComment 2) nor our patient stratification based on the proposed signatures (Extended Data Figure 13) differ\nsignificantly across treatment groups.\n15\na, Kaplan-Meier plots showing the association between PRRX1 expression and RFS using data sets from the Chin\nbreast adenocarcinoma series (Chin et al., 2006) and from a series of lung SCC (Raponi et al., 2006). b, Statistical\nassociation between PRRX1 with metastasis in (Chin et al., 2006). c, Representative examples of ISH for PRRX1 in\na retrospective grade 3 infiltrating ductal breast carcinoma series (n = 113) included in a tumor microarray\n(TMA). Lower panels show high power images. Scale bars: 200 μm (upper panels) and 50 μm (lower panels).\nStatistical association between PRRX1 expression and metastasis formation from the expression analysis of\nthe TMA. ∗n (%), number of cases (percentage); ns, nonsignificant.\nExtracted and modified from: https://doi.org/10.1016/].ccr.2012.10.012.\n9) Regarding hermetic role of Prrx1, is there a possibility that it may be modified by treatment?\nWe thank the reviewer for raising the important question of whether the hormetic role of Prrx1\nmight be\nmodified by therapeutic interventions. Conceptually, this is indeed plausible. Many standard treatments,\nsuch as cytotoxic chemotherapy, radiotherapy, or CDK4/6 inhibitors, directly influence cellular\nproliferation, which constitutes one of the key axes through which intermediate Prrx1 levels enhance\n16\nmetastatic competence. In principle, an antiproliferative therapy could attenuate or even eliminate the\nselective advantage of the Prrx1-intermediate, highly proliferative population relative to the Prrx1-high, less\nproliferative group, thereby altering the hormetic relationship with metastatic potential we describe.\nA related question is whether Prrx1 heterozygous tumors (Prrx1+/f) may exhibit differential sensitivity to\nantiproliferative treatments compared with Prrx1+/+ tumors, particularly regarding metastatic potential.\nAddressing this rigorously would require extensive therapeutic experiments in vivo, which would constitute\nan independent research project and is beyond the scope of the present study. We fully agree that this\nrepresents a highly relevant future direction.\nNevertheless, given the clinical interest of this point, we explored whether any treatment-associated\nmodulation of PRRX1 levels could be inferred from the METABRIC dataset. As shown in our previous\nanswer to comment 2, we stratified patients according to receipt of chemotherapy, radiotherapy, or\nhormone therapy, and compared PRRX1 expression between treated and untreated groups for each\nmodality. We did not observe consistent or significant differences in PRRX1 levels between these groups.\nHowever, because most patients receive at least one form of treatment, and untreated patients are\nextremely rare, it is difficult to isolate the effect of any individual therapy. If a particular treatment elevates\nPRRX1 levels across patients, such an effect would not be detectable simply by removing a single\ntherapeutic category from the comparison.\nWith a related aim, we examined whether our proposed combined invasion–proliferation stratification\n(based on biological features regulated by Prrx1) is preserved across treatment subgroups. These\nanalyses have been incorporated in the revised version of the manuscript as Extended Data 13 and\ndiscussed in page 12 of main text. Interestingly, the prognostic hierarchy of these groups was highly\nconserved across all treatment categories. In every treatment group, the Invasion-high / Proliferation-\nhigh group consistently associated with poor survival, whereas the Invasion-low / Proliferation-\nhigh group demonstrated a favorable prognosis. These findings suggest that our stratification strategy is\nrobust and broadly generalizable, and that its prognostic value is not driven by a particular therapy-\nrelated subgroup.\nOverall, while treatment may theoretically influence the hormetic behavior of Prrx1, our analyses indicate that\nat the patient-data level, PRRX1 expression and Prrx1-associated biological stratification appear remarkably\nstable across treatment contexts. Further dedicated experimental work will be essential to dissect treatment–\nPrrx1 interactions mechanistically.\n10) How does the addition of cell cycle regulators, including CDK4/6 inhibitors, to the BT549 cell line\n(with or without hormone therapy) affect these levels? This will likely be an important point in this study.\nWe agree with the referee on the interest of exploring the impact of treatments on Prrx1 levels and\nassociated signatures. Our preliminary analyses (as discussed in previous responses) indicate that Prrx1\nexpression and its related signatures are not significantly affected by standard therapies in the datasets\nanalysed.\nRigorous experimental evaluation of the effects of cell cycle regulators, such as CDK4/6 inhibitors, or hormone\ntherapies would require extensive in vitro and in vivo studies, as the hormetic role of Prrx1 can\n17\nonly be fully assessed within the context of the complete metastatic cascade. This constitutes a separate\nresearch project, which we may consider to pursue in the future. Nonetheless, based on our current\nunderstanding, we predict that CDK4/6 inhibition could partially reduce the enhanced metastatic potential\nobserved in BT549 cells following Prrx1 depletion. Regarding hormone therapy, BT549 cells are derived\nfrom a triple-negative breast cancer (TNBC) patient, lacking estrogen receptor, progesterone receptor,\nand HER2 expression, and therefore are not an appropriate model to study hormone-based effects.\n11) Is it possible to directly confirm the regulatory role of PRRX1 through additional functional assays\n(e.g., knockdown or rescue experiments targeting dormancy-related genes such as GAS6/MME)? We\nthank the reviewer for this insightful question regarding the possibility of performing additional\nfunctional assays, such as knockdown or rescue of dormancy-related genes (GAS6, MME), to further\nvalidate the regulatory role of PRRX1.\nThe identification of true functional regulators of dormancy is indeed a highly active and exciting field.\nHowever, regarding the specific dormancy-related genes included in the signatures used in our study (Ren\net al. Breast Cancer Research 2022 10.1186/s13058-022-01503-5), the authors of the original dormancy\nmanuscripts have already performed extensive gain- and loss-of-function experiments. In those studies,\noverexpression of dormancy-associated genes in non-dormant cell lines, or knockout in dormant cell lines\n(including GAS6 and MME), did not alter metastatic potential. These results indicate that these genes act\nas markers of a dormant transcriptional state, rather than as upstream drivers or effectors capable of\nmodulating dormancy on their own.\nDormancy is fundamentally defined by cell-cycle arrest, and accumulating evidence points to cell-cycle\nregulators as the primary functional controllers of this state. In line with this, our work identifies a specific\ncancer-cell subpopulation in vivo in which Prrx1 levels tightly associate with cell-cycle inhibition,\ndormancy-related transcriptional programs, and metastatic competence. Moreover, we demonstrate that\nmanipulating Prrx1 expression directly modulates these programs, including the activation of CDKN2\ninhibitor genes, which constitute well-established effectors of dormancy.\nTaken together, our data support the conclusion that Prrx1 itself functions as a bona fide regulator of\ndormancy, rather than acting downstream of previously described marker genes such as GAS6 or MME. For\nthis reason, additional manipulation of individual dormancy-signature genes is unlikely to provide further\nmechanistic insight into the core regulatory relationship between Prrx1, cell-cycle arrest, and metastatic\nbehavior.\nAs an additional functional assay, we have developed a new model by engineering MDCK cells to express\ngraded levels of Prrx1 and performed invasion and growth assays, as well as RT–qPCR analysis of cell-cycle\nregulators and dormancy-associated genes (Figure 6d-g and page 11 of main text).\n12) Intracellular metabolism, particularly cellular energy metabolism, appears to be important in the\nrelationship between Prrx1 expression levels and EMT, invasion, proliferation, and dormancy. If authors have\nany knowledge on this point, it will strengthen the hypothesis that the authors address.\nWe thank the referee for highlighting the potential link between intracellular metabolism and Prrx1-mediated\nphenotypes. Indeed, there is a well-established connection between EMT/invasion and\n18\ncellular metabolism, particularly glucose metabolism (Youssef and Nieto., Dev Cell 2020\n10.1016/j.devcel.2020.03.021). However, we did not present any data in this manuscript addressing the\nrelationship between Prrx1 expression levels and metabolic regulation, and it would be speculative to draw\nconclusions in this regard. However, we acknowledge it is an interesting avenue for future investigation.\n13) In addition, how have the authors considered the interactions between tumour cells with different Prrx1\nexpression levels at the invasion frontier? For example, what about energy exchange? What about interactions\nbetween tumour cells and stromal cells?\nWe agree with the referee that cell–cell interactions are important drivers of cellular behavior, particularly\nat the invasive front where cancer cells are intermingled with stromal populations. However, in this study,\nwe have focused on the cell-autonomous functions, as Prrx1 is a transcription factor, and we are specifically\naddressing the impact of different Prrx1 levels in cancer cells on metastatic potential.\nMoreover, due to technical limitations and the very low expression levels of Prrx1, it is currently not\nfeasible to reliably distinguish intermediate from high Prrx1–expressing cells in tissue sections. We\nacknowledge that future advances in spatial transcriptomics/metabolomics could allow a more detailed\nanalysis of these interactions, including metabolite exchange and tumor–stromal crosstalk.\n14) From the perspective of Prrx1 expression inducers or inhibitors, what conditions would result in\ndormancy break? It will be better to explain, even as a hypothesis. The relationship between various\ntreatments, such as endocrine therapy and anti-HER2 therapy, and dormancy break has been shown\nto be a long-term clinical phenomenon, in meta-analyses that have compiled many prospective cases.\nThis may lead to the development of new treatments that take into account the regulation of PXRR1\nexpression.\nWe appreciate the opportunity to discuss the potential translational implications of our findings regarding\nPrrx1 regulation and dormancy. Although speculative, our results suggest a model in which Prrx1 inhibitors\ncould promote the break of dormancy in already disseminated cancer cells, whereas Prrx1 inducers could\nfavor the maintenance of dormancy, potentially restraining metastatic outgrowth.\nClinically, long-term dormancy break is particularly relevant in hormone receptor–positive breast\ncancers, whereas it is less common in hormone receptor–negative subtypes. Some studies further\nsuggest that anti-HER2 therapies may promote dormancy in subsets of disseminated cancer cells,\nwhich could serve as seeds for long-term recurrence (Ruth et al., Breast Cancer Research 2021\n10.1186/s13058-021-01416-9. It is important to note, however, that Prrx1 is a transcription factor\nexpressed at very low levels and in multiple cell types beyond cancer cells. These characteristics\ncurrently limit its direct druggability. Nevertheless, emerging approaches, such as PROTACs or\nDUBTACs, could open avenues to target transcription factors like Prrx1 in the future. Indeed, we are\ncurrently writing a review with a section dedicated to anti-EMT future therapy development where we\ndiscuss the above-mentioned strategies among others.\n15) Statistical Robustness: As indicated, in some figures (e.g., metastasis incidence by PRRX1 level), the\nsample size is small. The authors should provide confidence intervals etc.\nWe thank the reviewer for the suggestion to improve the statistical robustness of our data. In response,\nwe\nhave included a new Supplementary Table 4 including pairwise comparisons between the patient\n19\ngroups shown in Figure 1. For each comparison, we report the odds ratio (OR) and relative risk (RR)\nalong with the corresponding 95% confidence intervals and p-values. Supplementary Table 4 also\nincludes similar analysis regarding stratification of breast cancer patients in Figure 7b.\nMinor Concerns\n0) Regarding the nonlinear correlation between PRRX1 expression level and metastasis incidence in\nbreast cancer patients, the paper states that, while TNBC is known to include multiple subtypes, this\ndata suggests that the high PRRX1 expression group may specifically belong to a TNBC subtype with a\nlow rate of metastasis. However, in Figure 1C, it is unclear how many TNBC cases actually existed among\nthe basal cases with high PRRX1 expression (n=8) and which TNBC molecular subtype they belonged\nto.\nWe thank the reviewer for this insightful comment, which allowed us to identify an inconsistency in the\nterminology used in the manuscript. In the text, we referred to TNBC, whereas the actual classification\nwas Basal. We have corrected this terminology throughout the manuscript.\nThe statement that high PRRX1-expressing Basal patients may belong to a “good prognosis TNBC\nsubtype” was intended as a hypothesis rather than a firm claim. We have included in the new version of\nthe manuscript new analysis of the METABRIC dataset revealing that the highest PRRX1 expression is\nobserved in the Claudin-low subtype (Extended Data Figure 1d-f and discussed in page 4 of main text) \nconsistent with its classification as a high-EMT subtype. Within Claudin-low patients, the majority belong\nto Integrative Cluster 4, which we identified as a genomic subgroup with highest PRRX1 expression\n(Extended Data Figure 1c) and it is considered a good prognosis subgroup (Curtis et al , Nature 2012\n10.1038/nature10983). In addition, recent detailed analyses of Claudin-low patients have revealed\ninternal heterogeneity, including a subgroup with very high EMT levels and low proliferation with\nparticularly good prognosis, which aligns closely with our findings (Pommier et al., Nat Comm 2020\n10.1038/s41467-020-17249-7) (see Figure next page).\nWhile this is beyond the fundamental biology focus of the present work, further detailed analysis of our\npatient cohort could clarify whether high PRRX1-expressing patients correspond to this good prognosis\nClaudin-low subgroup.\n1) In the METABRIC database analysis, the negative correlation with growth appears weak.\nHow definitive can this be?\nWe agree with the reviewer that the negative correlation between PRRX1 expression and proliferation\nscores in the METABRIC bulk RNA-seq dataset appears modest (Figure 7a), especially when compared\nwith the much stronger correlations observed with invasion or dormancy. However, we believe this weak\nassociation is entirely consistent with the biology we describe in this study:\n• First, our mouse models show an analogous pattern: Prrx1 deletion does not alter overall tumor\nburden or global tumor proliferation (Main Figure 2b-c). The impact on proliferation emerges only\nwhen focusing specifically on the invasive cancer-cell population, as demonstrated both by scRNA-\nseq and by immunofluorescence in invasive regions (Main Figure 4g-i). This is expected, as Prrx1\nis only expressed in invasive cells and accounts for <1% of the tumor cells (Figure 4e). Thus, a\nstrong negative correlation in bulk analyses would not be expected, as the proliferative phenotype\nof the vast majority of Prrx1 negative cells would dominate the signal.\n20\na, ssGSEA score distribution for each molecular subtype (including claudin-low subtypes), from METABRIC\nand TCGA cohorts, related to EMT. b, ssGSEA score distribution for each molecular subtype ((including\nclaudin-low subtypes) from the METABRIC cohort related to proliferation and estrogen response. c, Overall\nsurvival and disease-specific survival each molecular subtype (including claudin-low subtypes) in METABRIC.\nExtracted and modified from: https://doi.org/10.1038/s41467-020-17249-7\n• We believe a similar dilution effect occurs in the METABRIC dataset: bulk RNA-seq captures\ntranscriptional averages across thousands of heterogeneous tumor cells, making subtle\nproliferation-associated features of a minor Prrx1+ subpopulation difficult to detect at the whole-\ntumor level.\n• Additional evidence supporting this interpretation comes from the metastatic lesions in our mouse\nmodel. Although metastases from Prrx1 +/+ and Prrx1+/f tumors do show significant differences\nin proliferation (Figure 4l), these differences are clearly smaller than those observed in the\ninvasive regions of primary tumors (Figure 4i). This is likely due to required downregulation of\nPrrx1 expression for metastatic outgrowth (Ocaña et al., Cancer Cell 2012\n10.1016/j.ccr.2012.10.012), which again dilutes the correlation between Prrx1 levels and\nproliferative behavior.\ndoes not contradict our findings; rather, it reflects the spatial and cellular specificity of Prrx1-associated\nproliferation changes, which become evident only when the appropriate cell subpopulations are examined.\nReviewer 5\nWe thank the reviewer for the comments and appreciate the suggestions for improvement in the presentation\nof some of the data. We were surprised to read that they found our study descriptive when we have generated\ndifferent loss of function Prrx1 models with graded activity both in vitro and in vivo. We hope and believe that\nthe additional models presented in this revised version, and our response to all the specific questions and\nsuggestions including the additional quality controls are sufficient to overcome the concern.\n1) Fig1 and Fig2, the use of IHC provides important insights into protein expression; however, the absence\nof RNA-level validation (e.g., RNAscope) limits the ability to determine whether the observed protein changes\nare transcriptionally regulated. Incorporating RNA-based analyses, such as single-cell or spatial\ntranscriptomics, would significantly strengthen the conclusions.\nWe would respectfully like to clarify that our study does not deal with Prrx1 expression regulation, neither\nin patients nor in our mouse models but rather to teh regulation exertedby Prrx1 on its downstream targets.\nNevertheless, RNA-based analyses of Prrx1 expression have been reported by us previously (Ocaña et al.,\nCancer Cell 2012 10.1016/j.ccr.2012.10.012) (see Figure in our answer to comment 8 of reviewers 3/4; page\n16 of this response letter), where we were able to distinguish patient groups based on Prrx1 RNA scores.\nHowever, this approach didn’t allow to clearly distinguish the source of the signal (stromal or cancer cell).\nIn addition, we now include an extended plot in Figure 7d of the manuscript showing clear differences in\nPrrx1 transcripts levels across the different METABRIC groups.\nRegarding Figure 2, although we use IF to identify Prrx1+ cancer cells within invasive regions, we\ndemonstrated in previous work that deletion of Prrx1 alleles leads to a direct and measurable reduction\nof Prrx1 in embryonic derived fibroblasts at protein level (see second Figure in our answer to Reviewer\n2- Comment 2; page 6 of this response letter) and here we prove it at the RNA level in FACS-isolated\ncancer cells (Figure 2a) in the different allelic situations. We have also provided scRNA-seq data showing\nthat Prrx1 gene expression is restricted to a specific cancer cell population (Figure 4e). This same\npopulation is identified as invasive through spatial transcriptomics (Figures 2 and 4b), another approach\nmentioned by the reviewer, and corresponds to the regions where Prrx1 protein is detected (Figure 2f).\nAltogether, these data suggest that while post translational mechanism could play role, there seems to\nbe a strong correlation between the available alleles, and RNA and protein levels in the case of Prrx1.\n2) Extended Fig1, the authors should show each individual replicate for the cell boundary staining with\nDAPI, including zoomed-in views to demonstrate clear segregation of cells. In addition, quality statistics\nsuch as average detected genes (nGene) and counts (nCounts) for each replicate should be provided. We\nagree with the reviewer on the importance of providing a detailed assessment of data quality, particularly\nfor spatial transcriptomics, where accurate cell segmentation is essential. Following these\n22\nsuggestions, we have included the individual replicates for the cell boundary staining with DAPI in the\nExtended Data Fig. 3b, with high-magnification views to clearly illustrate cell separation. Additionally, we now\nprovide quality control metrics, including average detected genes (nGene) and total counts (nCounts) for each\nreplicate, in Extended Data Fig. 3a.\n3) ATAC-seq quality controls, Key metrics are missing, including unique fragments per cell, fraction of reads\nin peaks (FRiP), transcription start site (TSS) enrichment, fragment size distributions, and mitochondrial gene\nproportions. These are essential for assessing dataset quality.\nFollowing the reviewer’s comment, we have generated a new Extended Data Fig. 8 that includes the requested\nscATAC-seq quality control metrics. Thanks for the suggestion.\n4) Fig 5a-f, when plotting the marker gene peaks, all populations should be visualized to show its specificity\nof the population of interest.\nWe appreciate the referee’s effort to strengthen our scATAC analysis. We understand that the comment\nrefers to the plots showing Prrx1-binding to accessible regions in genes of interest (Figure 5e,f). We\napologize, as we did not mean that these accessible regions are specific to cells in cluster 15, but rather\nthat we focused our chromatin analysis on cluster 15 because it represents the Prrx1+ cancer cell\npopulation. We agree that the original text was not clear, and is now revised (page 9). Please see figure\nin next page showing the requested plots displaying chromatin accessibility across all cancer cell\npopulations. Except for Ccnd2, the identified regulatory regions are not uniquely accessible in cluster 15.\nHowever, this is not inconsistent with our interpretation. These results indicate that the specificity of\nPrrx1-mediated regulation does not arise from chromatin accessibility in specific cells clusters, but rather\nfrom the restricted expression of Prrx1 in cluster 15 cells.\nPeaks with identified binding sites for Prrx1 are indicated in red.\n5) Fig 5i, the ChIP-seq IGV tracks of different cell lines should be normalized and visualized on the same\nscale to allow direct comparison.\nWe fully agree with the referee that, for direct comparison, ChIP-seq tracks from different cell lines\nshould ideally be visualized on the same scale. However, it is important to clarify that the objective of this\nanalysis was not to compare signal intensities quantitatively across cell lines, but rather to determine whether\nbinding peaks could be identified in breast cancer cell lines expressing Prrx1. Indeed, these cells have different\nPRRX1 levels. In this context, using an adjusted scale for each sample allows for optimal visualization of peaks\nrelative to background.\nNonetheless, in response to the referee’s suggestion, we have generated here a version of the IGV tracks\nnormalized to the same scale across all cell lines (see Figure below). We would be pleased to include it as\nExtended Data or to replace the current version if the reviewer deems it necessary.\ninvasion genes in 3 human breast cancer cell lines: HS578 in blue, MDAMB157 in green, MDAMB436 in red.\n6) Detailed Code should be provided.\nThanks for the suggestion. We have prepared and uploaded the complete and detailed code pipeline to\na Zenodo repository (see link below). The repository is currently accessible to the reviewers through\nthis private link and will be made publicly available upon publication.\nhttps://zenodo.org/records/17726450?token=eyJhbGciOiJIUzUxMiIsImlhdCI6MTc2NDI0MjY4NSwiZX\nhwIjoxNzc3NTA3MTk5fQ.eyJpZCI6IjBlNjQyOTRhLTVmMGItNDU5NC05YjFhLTM3ZjY4YjI4ODEwMiIsI\nmRhdGEiOnt9LCJyYW5kb20iOiI2NDQ0MmZmY2ZkYzFkNDdjNzAyZjBjMDljZWY3N2U4YiJ9.nbnL06LJ\nrckNLBJ7Jb2_THfvejl7aVys0tx-gmJTNY7aFB6iVnfhaXf_7x2jDRO6shsEgKC7QHsq-DE0rVHl2g\nWe thank the reviewer for acknowledging the use of different experimental models and for the comments on\nthe novelty, quality and robustness of our data. Here we are pleased to provide additional details of data\nanalysis as requested.\n1) Fig3f - there are some stromal genes among the upregulated genes, such as Col1a1, Postn, Tnc - how\nsure are you that those are cancer cells and not CAFs?\nWe thank the reviewer for highlighting this important point, which is indeed critical when interpreting gene\nexpression changes in EMT-associated contexts. As noted by the reviewer, genes such as Col1a1, Postn,\nand Tnc, highly expressed by CAFs, are upregulated in cancer cells located within invasive areas when compared\nto cancer cells in non-invasive regions. This observation is in line with well-established EMT biology, as\nactivation of mesenchymal programs is a hallmark of invasive EMT states.\nImportantly, the differential expression analysis presented in Fig. 3f was performed exclusively on cancer\ncells, as defined by the clustering annotations derived from the full gene panel. Nonetheless, we fully\nagree that the close physical proximity of cancer cells and CAFs in invasive regions raises the possibility\nof transcript misassignment or miss-segmentation. To further validate that the expression is in cancer\ncells, we performed an additional inspection of representative invasive and non-invasive areas. Our spatial\ntranscriptomics custom panel also includes qualitative tdTomato probes (not used for clustering) that\nallow unambiguous discrimination between tumor and stromal compartments (Extended Data Fig. 2). We\nnow include an example of the analysis in a new figure (Extended Data Fig. 4), which confirms that:\n1. Stromal cells, especially CAFs, show the highest absolute expression of Col1a1, Postn, and Tnc, as\nexpected.\n2. Cancer cells identified both by clustering and tdTomato expression show clear upregulation of these\ngenes specifically within invasive areas, coinciding with high Krt14 expression.\n3. Cancer cells in non-invasive regions show little to no expression of these genes (not included in\nExtended Data Figure 4, see Figure below).\n4. Consistently, in our scRNA-seq data, the Prrx1+ cancer cell cluster expresses all three genes at\nhigher levels than other tumor cell clusters.\nTogether, these analyses confirm that the genes mentioned by the reviewer are genuinely upregulated in\ninvasive EMT-high cancer cells reflecting the activation of mesenchymal and ECM-remodeling programs\ncharacteristic of EMT and invasion.\nexpression. The images show the expression of mesenchymal markers (dots) and tdTomato (green),\nas in\nExtended Data Fig. 4.\n2) Fig 4 - Cluster 25 seems to be the one with the highest Prrx1 expression. Is this enriched in a specific\ngenotype, from extended figure 3c it doesn’t see so, but can you show one umap coloured by genotypes, or\n3 maps as in extended figure 3c but coloured by cluster to make it clearer?\nWe thank the reviewer for this suggestion, which indeed helps clarify the relationship between Prrx1-high\nclusters and genotype distribution. As requested, we now show UMAPs coloured by cluster (Extended Data\nFigure 6c).\n3) Fig 4 - Cluster 20 seems to have almost no Prrx1 expression, but still has high EMT and invasive score,\nhow do you explain this?\nWe thank the reviewer for this insightful comment. We agree that cluster 20 displays high EMT and invasive\nscores despite showing almost no detectable Prrx1 expression, and this is indeed an observation we have\ncarefully considered. Firstly, we acknowledge the possibility that Prrx1, being a very lowly expressed\ntranscription factor, may be present in cluster 20 at levels below the detection limit of scRNA-seq.\nNevertheless, we consider it unlikely that Prrx1 plays a major functional role in this cluster, as we do not\nobserve differences in downstream transcriptional programs (epithelial, mesenchymal, invasion, proliferation;\ndata not shown) across genotypes within cluster 20. This contrasts with cluster 25, where genotype-\ndependent differences in these signatures are clearly detectable.\nhave not reached a critical threshold to translate it into effective invasive behavior. In a recently\npublished work in in the same mouse model (Youssef et al., Nat Cancer 2024 10.1038/s43018-024-00839-5),\nwe showed that these clusters lie along a continuous EMT trajectory, and other EMT\ntranscription factors such as Zeb1 and Twist family members are also expressed along this trajectory.\nThese TFs may be sufficient to initiate EMT and invasion-related programsbut, in the absence of Prrx1,\nthey seem to fail to drive the full functional activation characteristic of cluster 25.\n4) For the scRNA-seq analyses, the authors write that “The top 4000 HVGs were selected for the Prrx1\nf/f and Prrx1 +/f samples, whereas top 6000 HVGs were selected for Prrx1 +/+ samples”. Can they\nexpand on why this was done?\nWe thank the reviewer for the opportunity to clarify this point. The choice of selecting 6000 HVGs for the\nPrrx1 +/+ samples and 4000 HVGs for the Prrx1 +/f and Prrx1 f/f samples was motivated by technical\ndifferences in sequencing chemistry across datasets.\nThe Prrx1 +/+ samples were generated using a combination of 10x Genomics v2 and v3 chemistries,\nwhereas the Prrx1 +/f and Prrx1 f/f samples were sequenced exclusively with v3. As documented by 10x\nGenomics (Technical Note, 2019), v3 chemistry yields substantially higher mRNA capture efficiency and\ndetects a greater number of genes per cell compared to v2 (see Figure below). To mitigate this\ndiscrepancy in gene detection sensitivity and to ensure a comparable representation of biological\nvariability, we calculated HVGs per sample and then subsequently took their intersection by condition.\nThis strategy allowed us to balance the number of informative genes across conditions while limiting the\nintroduction of sample specific technical noise.\nBoxplots showing number of detected genes (left) and features (right) across Prrx1 +/+ samples coloured by\nchemistry used to prepare the library.\nBecause v2 samples\ndetect fewer genes\noverall, we increased\nthe number of HVGs\ncomputed for the\nPrrx1 +/+ dataset\n(6000 HVGs) to avoid\nlosing informative\ngenes that are\nrobustly captured in\nv3 samples of the\nsame condition (as we are only taking the commonly detected in all the samples from the same condition).\nFor the Prrx1 +/f and Prrx1 f/f samples, which were uniformly processed with v3, we selected\n28\nspace, with a very similar number of HVG coming from each genotype (which would not happen if 4000\nHVGs were calculated for Prrx1 +/+ condition). This strategy was designed together with experts in atlas\ngeneration and data integration (Holger Heyn´s group at CRG, Barcelona) as indicated in\nacknowledgements.\nImportantly, as detailed in the Methods, although all samples contributed to clustering, all differential\ngene expression analyses were restricted to samples sequenced with the same chemistry (v3): Prrx1\n+/+ (3 samples), Prrx1 +/f (3 samples), and Prrx1 f/f (4 samples). This ensured that genotype\ncomparisons were not confounded by differences in sequencing chemistry.\n5) Fig 5 - the authors imputed gene expression in scATAC-seq data using the scRNA-seq as a\nreference. Do they see a similar enrichment of the STINS score if they use the gene activity\nscores instead of imputation?\nWe thank the reviewer for this constructive suggestion. In response, we have performed the requested\nanalysis by calculating the STINS score directly from gene activity values, rather than from imputed gene\nexpression derived from scRNA-seq. As shown in the Figure below, although the dynamic range of the\nsignal is reduced, an expected outcome given the lower resolution of gene activity scores, we still observe\na clear enrichment of the STINS signature in clusters 14 and 15, matching the pattern obtained with the\nimputed expression approach.\nThis consistency supports the robustness of our findings and indicates that the chromatin regions\nassociated with promotors of STINS-related genes are indeed more accessible in these specific cell\nstates, reinforcing the conclusions drawn from the original analysis.\na, Enrichment of STINS (spatial transcriptomic-derived invasive signature) gene activity score represented over\nscATAC UMAP. b, Boxplot of STINS gene activity score across cancer cell clusters in scATAC data.\nMinor:\n6) As someone who\ndoesn't work with\nbreast cancer, it was\nnot clear why the\nauthors used Krt14\nareas at the tumour\nedge to identify\ninvasive areas, perhaps a very brief explanation could be added.\n29\n• We apologize for not making that clearer. In the PyMT model, we and others have previously\nshown that, during tumor progression, luminal tumor cells undergo a basal-like transition\ncharacterized by de novo Krt14 expression, which coincides with the activation of an invasive\nEMT program (Youssef et al., Nat Cancer 2024 10.1038/s43018-024-00839-5). This is\nconsistent with additional studies demonstrating that Krt14 marks the invasive and\ndisseminative population in this cancer model (Cheung et al, Cell 2013\n10.1016/j.cell.2013.11.029).\nBoth in this manuscript and in our prior work characterizing the invasive EMT program in detail (Youssef\net\nal., Nat Cancer 2024) Krt14 expression was used to identify invasive areas in combination with the presence\nof disrupted or unstructured laminin borders, a well-established morphological hallmark of invasive behavior.\nWe have clarified this in the revised manuscript (see page 5).\nReviewer #3/4 (Remarks to the Author):\nThe authors have responded to the series of comments with commendable thoroughness. They have\nundertaken substantial additional analyses, experiments, and reflections, and these efforts not only greatly\nfacilitate the understanding of the present findings but also deepen the overall significance of the study.\nThe rationale for restricting the analysis to grade 3 TNBC has now been supported by careful\nargumentation and additional data.\nFinally, I would ask the authors to explicitly summarise, within the main text, the key directions for\nfuture research that emerge from this work.\nWe thank the reviewer for acknowledging our efforts to improve the manuscript. As discussed during\nthe revision process, we agree that including additional histological subtypes and a larger number of\npatients will be important to enhance the translational relevance and generalizability of these findings.\nAccordingly, we have added a new sentence in the Discussion (page 14) mentioning key directions for\nfuture research.\nReviewer #5 (Remarks to the Author):\nIn the revised version of the manuscript, M. Angela Nieto and colleagues addressed some of my\nconcerns in detail with convincing data and explanation. While the mechanistic link between PRRX1\nand its targets remains to be strengthened:\n- Line303-305, the authors mentioned that “Transcription factor binding motif analysis showed that\nPrrx1 binding sites are more accessible in cluster 15 compared to the rest of the clusters (Fig. 5d),\n(P<1.93*10^-18)” However, they only found Ccnd2 is uniquely accessible to cluster 15. Given the\nsignificance of this motif accessibility in cluster15, did they look at what other target genes unique to\ncluster15 and associated functions of these target genes?\nWe appreciate this comment. This is clarification that does not require any modification of the text.\nAs\nthe reviewer notes, the transcription factor binding motif analysis revealed that Prrx1 binding sites are\nsignificantly more accessible in cluster 15 compared to the other clusters (Fig. 5d). This analysis was\nperformed genome-wide and was intended to confirm the expected regulatory role of Prrx1 as a\ntranscription factor, as well as its higher specificity in cluster 15, where Prrx1 is predicted to be\nexpressed.\nImportantly, Ccnd2 is not the only gene showing increased accessibility in cluster 15. Rather, Ccnd2\nis the only gene among the six selected candidates shown in Fig. 5e that displays higher accessibility\nin this cluster. In fact, multiple genes harboring Prrx1 binding sites exhibit increased accessibility in\ncluster 15, particularly those related to EMT, the canonical cell program regulated by Prrx1. Critically,\nwe did not claim in the manuscript that the selected genes are more accessible in this cluster compared\nto rest.\nIn addition, the goal of the analysis presented in Fig. 5e was not to provide a comprehensive\ncharacterization of all Prrx1 target genes. At that point in our study, after the scRNA-seq analysis (Fig.\n4), the scope had already been narrowed to a focused set of genes related to cell cycle regulation and\ndormancy.\n- Page9, line 312: The authors claimed that Prrx1 can plays an important regulatory role, as it can\ndirectly bind promoters of genes associated with cell-cycle control and dormancy. However, the\nmechanistic links are not sufficiently established here, especially when they want to claim direct\nregulation of dormancy (not unique to cluster 15).\nThis is clarification that does not require any modification of the manuscript.\nRegarding dormancy, while it is not strictly unique to cluster 15 in the scATAC-seq data, our scRNA-\nseq data show that the dormancy signature is highly specific to cluster 25, the only Prrx1-positive\ncluster (Fig. 4e), and is reduced upon Prrx1 downregulation in Prrx1 +/f tumors (Fig. 4k). Members of\nthe CDKN2 family show a similar enrichment pattern (Fig. 4]). Thus, in the case of cell-cycle and\ndormancy genes, regulatory specificity appears to derive from Prrx1 expression rather than differential\nmotif accessibility.\nWe find open chromatin regions enriched for Prrx1 binding motifs in the promoters of target genes\nwithin the only cluster where Prrx1 is expressed.\nThe author may need to perform gold standard assays such ChIP-seq or CUT&RUN (Peter J., et al,\nelife, 2017) to make this claim more solid.\nThis is the point you mentioned that the requested experiments lie out of the scope of this work. We\nagree and thank you very much. While we acknowledge that ChIP-seq and CUT&RUN are gold-\nstandard methods to demonstrate direct binding, their application in this in vivo context is technically\nextremely challenging. Prrx1-positive cells represent ≤1% of tumor cells, cannot currently be efficiently\nisolated,\nand express Prrx1 at very low levels, limiting signal detection. Moreover, these assays demonstrate\nbinding but not regulation. For these reasons, we used scATAC-seq to interrogate the regulatory\nlandscape in vivo by clustering Prrx1-positive cells, an approach not feasible by sorting.\nIn addition, the conclusion is also supported by multiple complementary approaches: Prrx1\ndownregulation and overexpression modulate cell-cycle and dormancy genes in the predicted direction\n(Fig. 6b, 6g), and ChIP-seq data from breast cancer cell lines demonstrate direct Prrx1 binding to\nseveral of these genes (Fig. 6c). Together, these data support a model in which Prrx1 directly regulates\nthe proposed targets.\n- In Extended Data Figure 8, the authors presented the unique fragments per cell. However, the median\nseems too good to be unique. Probably they may double check whether duplicates were removed? In\naddition, the FRiP is supposed to fall into 0-1. They may check (De Rop, Florian V., et al. Nature\nBiotech, 2024.) for these analyses.\nWe thank the reviewer for carefully assessing the quality of our scATAC-seq data. The values for\nunique fragments per cell are taken directly from the Cell Ranger output for 10x Genomics scATAC-\nseq, which includes duplicate removal. Although these values are higher than commonly reported, they\ndo not arise from a technical artifact. On average, we obtained ~44k raw reads per cell, with >95%\nhigh-quality fragments across all samples. Importantly, the data were generated using the 10x\nGenomics V2 kit, which, as shown in the reference cited by the reviewer (De Rop et al., Nature\nBiotechnology, 2024), yields higher numbers of unique fragments in peaks compared to other\nchemistries. This is clarification that does not require any modification of the manuscript.\nIn addition, sample preparation was extensively optimized to obtain highly pure and high-quality nuclei.\nThe data also derive from a spontaneous cancer mouse model with high cellular plasticity, reflected in\nthe large diversity of cell states observed in our scRNA-seq and scATAC-seq analyses. This plasticity\nmay be associated with globally increased chromatin accessibility, leading to a higher number of Tn5\ninsertion sites and fragments, and may also explain the moderate TSS enrichment, as accessibility\nextends beyond promoter regions. This is clarification that does not require any modification of the\nmanuscript.\nRegarding FRiP, we acknowledge an error in its initial calculation due to technical limitations of the\nSignac built-in function, which does not seem to perform correctly with 10x-derived data. We have now\nrecalculated FRiP directly from the Cell Ranger output and updated the plots accordingly. The revised\nvalues fall within the expected 0–1 range and are comparable to those reported for the V2 kit in (De\nRop, Florian V. et al. Nature Biotech, 2024) We thank the reviewer for identifying this issue, which has\nbeen solved.\n- Regarding visualization of the ChIP-seq IGV tracks, while using the same scale after sequencing-\ndepth normalization (e.g., CPM/RPM) enables fair comparison, interpretation remains highly\ndependent on antibody quality. In particular, insufficient discrimination between background and true\nsignal may confound the assessment of TF binding.\nWe also acknowledge this comment, that we address for you here and that does not require any\nchange in the manuscript. We agree with the reviewer that antibody quality is a critical factor for the\ninterpretation of ChIP-seq and CUT&RUN data, particularly when assessing transcription factor binding\nfrom IGV tracks. This concern was indeed one of the reasons why we did not rely on these approaches\nto study Prrx1 binding in vivo, where the analysis is substantially more challenging.\nWhile ChIP-seq data can be noisy, the datasets shown here display clear enrichment over background at\nseveral loci in multiple cell lines, allowing discrimination between signal and background in these cases.\nWe note that Prrx1 expression levels vary across cell lines, and therefore binding is not expected to be\nuniform across all targets or models.\nAs stated in the manuscript, the ChIP-seq datasets used in this study were previously published\n(Jovanovic et al., Cell Reports, 2023), and subjected to extensive quality control and integrative\nanalyses across multiple modalities. Importantly, in the present work these data are not interpreted in\nisolation but are integrated with our scRNA-seq and scATAC-seq data from in vivo tumors, genetic\nmodulation of Prrx1 in vivo and in vitro, and correlations in human datasets. Together, these\ncomplementary approaches support a regulatory role for Prrx1 despite the inherent limitations of ChIP-\nseq-based visualization protocols."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The reviewers have adequately addressed my concerns", "reviewer_2": "The authors have significantly improved the study and its impact; they also addressed my previous concerns.", "reviewer_3": "The authors have responded to the series of comments with commendable thoroughness. They have undertaken substantial\nadditional analyses, experiments, and reflections, and these efforts not only greatly facilitate the understanding of the present\nfindings but also deepen the overall significance of the study. The rationale for restricting the analysis to grade 3\nTNBC has now been supported by careful argumentation and additional data.\nWith regard to the limited sample size, I look forward to future studies conducted on a larger scale. The authors have also\nbegun to address, at least in part, the question of generalizability—namely, how these findings might apply to other\nsubgroups—which I find encouraging. Nevertheless, as with the previous point, I would strongly encourage further,\nmore\nin‑depth investigation of this issue in forthcoming work.\nHistological grade is undoubtedly important, but it is only one of several factors that inform therapeutic decision‑making.\nIn\nroutine practice, interobserver variability in grading among pathologists is a well‑recognized issue. Extending the analysis\nbeyond grade 3 TNBC to include other histological subgroups would, in my view, substantially strengthen the robustness\nand broader relevance of the present findings.\nFinally, I would ask the authors to explicitly summarise, within the main text, the key directions for future research that\nemerge from this work.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.", "reviewer_5": "In the revised version of the manuscript, M. Angela Nieto and colleagues addressed some of my concerns in detail with\nconvincing data and explanation. While the mechanistic link between PRRX1 and its targets remains to be strengthened:\n- Line303-305, the authors mentioned that “Transcription factor binding motif analysis showed that Prrx1 binding sites are\nmore accessible in cluster 15 compared to the rest of the clusters (Fig. 5d), (P<1.93*10^-18)” However, they only found\nCcnd2 is uniquely accessible to cluster 15. Given the significance of this motif accessibility in cluster15, did they look at\nwhat other target genes unique to cluster15 and associated functions of these target genes?\n- Page9, line 312: The authors claimed that Prrx1 can plays an important regulatory role, as it can directly bind promoters of\ngenes associated with cell-cycle control and dormancy. However, the mechanistic links are not sufficiently established here,\nespecially when they want to claim direct regulation of dormancy(not unique to cluster 15). The author may need to perform\ngold standard assays such ChIP-seq or CUT&RUN (Peter J., et al, elife, 2017) to make this claim more solid.\n- In Extended Data Figure 8, the authors presented the unique fragments per cell. However, the median seems too good to be\nunique. Probably they may double check whether duplicates were removed? In addition, the FRiP is supposed to fall into 0-1.\nThey may check (De Rop, Florian V., et al. Nature Biotech, 2024.) for these analyses.\n- Regarding visualization of the ChIP-seq IGV tracks, while using the same scale after sequencing-depth normalization (e.g.,\nCPM/RPM) enables fair comparison, interpretation remains highly dependent on antibody quality. In particular, insufficient\ndiscrimination between background and true signal may confound the assessment of TF binding.", "reviewer_6": "The authors have addressed all of my comments."}, "author_response": ""}]

Decision letter