Introduction

Tumor cell phenotypic heterogeneity is the product of genetic and epigenetic variability, localized differences in the microenvironment and stochastic fluctuations in intracellular signaling pathways (1). This diversity in cell phenotypes drives tumor evolution (2). Cells with traits that provide a competitive advantage over siblings become predominant, improving the fitness of the population (3). High variance in cell phenotypes also increases the odds that a subpopulation has the potential to survive when the local environment changes, such as during treatment or when disseminating cells seed other organs (4, 5). Interestingly, analysis of intrinsic and genetically engineered heterogeneity has revealed that phenotypic diversity also promotes functional interactions between tumor cells (6, 7). Rare enabler tumor subpopulations can promote the growth, immune evasion, invasion and metastasis of abundant siblings harboring different phenotypic traits (8-13). However, the regulatory programs that promote functional relationships between phenotypically distinct subpopulations remain poorly understood. This lack of mechanistic insight limits our ability to define the presence of functional tumor cell interactions, evaluate their impact on patient outcome, and develop new treatment options designed to target rare enabler cells (14).

To better understand how subpopulation relationships are induced during tumor progression, we sought to understand how tumor cells are conferred with the ability to promote the collective invasion of intrinsically less invasive siblings. Collective invasion is the predominant mode of invasion in carcinomas and promotes the dispersion of tumor clusters that seed metastases (15-19). Subpopulations of trailblazer cells detected in breast, lung and thyroid tumor populations have an enhanced ability to reorganize the extracellular matrix (ECM) into micro-tracks that facilitate the initiation of collective invasion (20-23). Importantly, sibling opportunist populations invade through the micro-tracks created by minority trailblazer populations in three-dimensional culture models and xenograft tumors (21, 24, 25). Moreover, opportunist cells can have greater proliferative capacity than trailblazers in vitro and in xenografts, indicating that trailblazer-opportunist relationships have the potential to increase the proliferative fitness of invading cell populations (24, 25). Thus, trailblazer cells promote phenotypic diversity in invading cell cohorts, which may increase the risk of metastasis and elevate the risk of resistance to treatment (26, 27).

Trailblazer cells are distinguished from the bulk tumor population by the increased expression of trailblazer genes, including regulators of the actin cytoskeleton, integrins, receptor tyrosine kinases, keratins, and cancer testes antigens (20, 21, 23, 24, 28). Trailblazer genes contribute to the formation of actin-rich cellular protrusions that provide traction and exert tensile forces necessary for ECM remodeling and other undefined functions (21, 23, 25). Interestingly, trailblazer genes are not required for other forms of cell movement, including migration along two-dimensional surfaces, movement within spheroids and opportunistic collective invasion (21). While genes that specifically contribute to unique invasive characteristics of trailblazer cells have begun to be defined, the regulatory programs underpinning the conversion to a trailblazer state remain largely unknown. Trailblazer cells in breast and lung cancer cell line models are epigenetically distinct from opportunist siblings, although the processes that initiate the conversion to the trailblazer state are not known (21, 23). Breast cancer trailblazer cells have basal characteristics, with Keratin 14 (Krt14) expression proposed as a universal feature for identifying trailblazer cells in mammary tumors (18, 20, 29, 30). Yet, the program controlling the induction of a basal trailblazer state is undefined. Breast cancer trailblazer cells also have features of epithelial-mesenchymal transition (EMT) re-programming (21). However not all EMT programs yield trailblazer activity and whether EMT programs regulate the expression of trailblazer genes is not known (25, 31). Thus, if and how the EMT process imbues a trailblazer state is unclear. Our poor understanding of how trailblazer populations are induced limits options for potential interventions targeting trailblazer cells and advances in employing subpopulation identification to refine patient outcome predictions.

To define factors regulating trailblazer program activity, we determined how cells were induced to collectively invade in the C3(1)/SV40 T-antigen (C3-TAg) genetically engineered mouse model (GEMM) of breast cancer. Based on histological analysis, breast cancer is defined as Estrogen Receptor-positive (ER-positive), Human Epidermal Growth Factor Receptor-positive (HER2-pos) or triple-negative (TNBC), which have low to undetectable levels of ER, HER2 and the Progesterone Receptor (32). Breast cancer is further stratified into luminal A, luminal B, HER2-enriched and basal-like based on mRNA expression. Luminal A are Estrogen Receptor-positive (ER-pos) while luminal B are more proliferative and are ER-pos or HER2-pos (33). The HER2-enriched subtype is generally HER2-pos and 70-80% of the basal-like breast cancers are classified as TNBC (33). The C3-TAg model is most similar to human basal-like TNBC (34, 35).

Through the analysis of tumor organoid invasion dynamics, a data driven functional dissection of cellular re-programming and quantitative imaging of fixed tumor sections, we discovered a Tgfβ regulated program that was necessary to induce and sustain a trailblazer state in C3-TAg tumors and validated this discovery in breast cancer patient tumor samples. Surprisingly, the Tgfβ regulated trailblazer cells were not specified by Krt14 expression and were more aggressive than a previously defined Krt14 expressing trailblazer population detected in a luminal B GEMM, indicating that there are multiple different programs for conferring trailblazer ability in breast tumors. In addition, trailblazer cells induced by Tgfβ signaling were themselves heterogeneous, having different degrees of re-programming and capacities to initiate invasion. Tgfβ signaling collaborated with a parallel Egfr-Mek1/2-Erk1/2-Fra1 pathway to coordinate the expression of a subset of pro-motility genes, including Axl. Interestingly, the integration was also equivalently antagonistic in scope, with Egfr activity restricting the Tgfβ induced expression of a cohort of pro-invasive genes. Notably, clinically approved inhibitors of Egfr, Mek1/2 and Axl suppressed collective invasion. This indicated that the collateral inhibition of subset of pro-invasive genes by the Egfr pathway was a feature of the trailblazer state. Thus, our results reveal how a new type of trailblazer cell is induced by a therapeutically addressable compromise between Tgfβ and Egfr signaling programs.

Results

The trailblazer cells in C3-TAg basal-like tumors are distinct from trailblazer cells in PyMT luminal B tumors

To define factors that initiate the conversion to a trailblazer state, we investigated the invasive properties of tumor cells isolated from the C3-TAg GEMM of basal-like TNBC (34, 35) and (Supplementary file 1). We first determined if there was heterogeneity with respect to the ability to invade in the tumor populations, which would indicate the potential for trailblazer-opportunist relationships. Indeed, there was extensive heterogeneity in the invasive ability of tumor organoids and single dissociated tumor cells (Figure 1A, Figure 1—figure supplement 1A-B, Figure 1—Videos 1-2). Clonal invasive (D6 and C7) and clonal non-invasive cells (B6 and G2) were also derived from a single C3-TAg tumor (Figure 1—figure supplement C-E, Supplementary file 1). To test whether the invasive subpopulations in C3-TAg tumors could lead the invasion of siblings, we generated heterogeneous clusters containing C3-TAg tumor cells and intrinsically non-invasive B6 cells (Figure 1B, Figure 1—Videos 3-5). The C3-TAg tumor cells and the D6 cells led the collective invasion of B6 cells (Figure 1B, Figure 1—figure supplement 1E, Figure 1—Videos 3 and 6). These combined results indicated that there were trailblazer subpopulations in C3-TAg tumors with an enhanced ability to lead the collective invasion of intrinsically less invasive siblings.

Figure 1 with 1 supplement.

To understand how the trailblazer state was induced in C3-TAg tumor cells, we next evaluated the relationship between Krt14 expression and the ability of cells to lead collective invasion. Krt14 identifies cells that have converted from a luminal non-invasive state to a basal trailblazer state in the MMTV-PyMT (PyMT) GEMM model of luminal B breast cancer, which is commonly used to study breast cancer invasion and metastasis (20, 36). Based on these results, it has been suggested that a basal gene expression program exemplified by Krt14 expression is a universal feature of breast cancer trailblazer cells (20, 29). Consistent with prior results, Krt14 was expressed in PyMT trailblazer cells and Krt14 expression correlated with PyMT tumor organoid invasion (Figure 1C-D). Krt14 was expressed in a greater proportion of cells in C3-TAg tumors compared to PyMT tumors (Figure 1C). The high Krt14 expression was consistent with the early acquisition of basal traits detected in C3-TAg and TNBC tumors during the initial stages of growth within mammary ducts (37, 38). However, Krt14 expression did not correlate with C3-TAg tumor collective invasion or a trailblazer state in our analysis of primary tumors, tumor organoids and trailblazer and opportunist cell lines (Figure 1C, E-F, Figure 1—figure supplement 1F). In fact, extensive collective invasion was detected in C3-TAg tumor organoids with few, if any, Krt14 expressing cells (Figure 1E). Krt14^neg C3-TAg tumor cells also led collective invasion as frequently as Krt14^pos cells (Figure 1F). C3-TAg organoids were relatively more invasive than PyMT organoids, demonstrating that the C3-TAg tumor invasion was more aggressive than the Krt14 specified invasion of PyMT tumors (Figure 1G, 1H, Figure 1—figure supplement 1G, Figure 1—Videos 7-8). Together, these results suggested that essential features of the C3-TAg trailblazer program were different than those of the PyMT trailblazer program, and that and more than a single Krt14 specified trailblazer program functioned in breast cancer patient tumors (Figure 1H).

Trailblazer cells in C3-TAg GEMM tumors are distinguished by a subset of EMT traits

Given the lack of concordance with the established Krt14 specified trailblazer program, we next took an unbiased approach to identify the key trailblazer regulatory program features in C3-TAg tumors. To do this, we used next generation RNA sequencing to define differentially expressed genes (DEGs) in tumor organoids from 4 C3-TAg tumors compared with non-invasive spheroids formed by 3 clonal C3-TAg tumor derived cell lines (Figure 2A, Figure 2—figure supplement 1A). Among epithelial keratins, Krt14 expression was reduced in C3-TAg tumor organoids compared to the non-invasive clonal C3-TAg cells, consistent with a lack of Krt14 involvement in C3-TAg collective invasion (Figure 2B, Supplementary file 2). Expression of the transcription factor p63 was also reduced in C3-TAg tumor organoids and trailblazer clones, further indicating that canonical traits of normal basal mammary epithelial cells and Krt14 specified trailblazer cells (20) were not essential features of C3-TAg trailblazer cells (Figure 2B, Figure 2—figure supplement 1B, Supplementary file 2). In addition, Krt14 expression was only detected in 1 out of 4 trailblazer populations in human breast cancer cell lines, which were all classified as TNBC (Figure 2—figure supplement 1C). The expression of p63 was also low in all 4 trailblazer populations (Figure 2—figure supplement 1C), consistent with our previous results showing that p63 expression is sufficient to suppress the induction of a trailblazer state (25).

Figure 2 with 1 supplement.

Genes highly expressed in C3-TAg tumor organoids were associated with biological processes that contribute to invasion (ECM Receptor Interaction, Focal Adhesions and Regulation of the Actin Cytoskeleton) and EMT (Figure 2C-D, Figure 2—figure supplement 1D-E, Supplementary file 2). The C3-TAg tumor organoids retained canonical epithelial gene expression as well, indicating that the majority of cells were in a hybrid EMT state (Figure 2E, Supplementary file 2). The EMT inducing transcription factors (EMT-TFs) Zeb1, Zeb2 and Fra1 were more highly expressed in C3-TAg organoids, consistent with their high expression in trailblazer cells in human breast cancer cell lines (Figure 2F, Figure 2—figure supplement 1C). However, the EMT-TFs Snail, Slug and Twist1 were expressed at similar or higher levels in the noninvasive clones, which was consistent with a lack of concordance between the induction of a trailblazer state and Snail, Slug and Twist1 expression in human breast cancer cell lines (Figure 2F, Figure 2—figure supplement 1C). These results indicated that the specific features of a trailblazer EMT program determined the invasive phenotype of C3-TAg tumor cells, rather than simply the presence or absence of EMT traits, consistent with results obtained analyzing human breast cancer cell lines (25, 31). Our combined results demonstrated that an EMT program exemplified by Vimentin expression conferred a trailblazer state in C3-TAg tumor subpopulations (Figure 2G).

To understand the relationship between the EMT program features of the C3-TAg organoids and collective invasion, we evaluated the expression of the microfilament protein Vimentin, a common EMT marker that was more highly expressed in C3-TAg organoids (Figure 2D). Vimentin expression was more frequently and extensively detected in invasive C3-TAg tumor cells compared to non-invasive C3-TAg mammary intraepithelial neoplasia (MIN) cells, consistent with Vimentin expression specifying a trailblazer EMT state in C3-TAg tumors. Notably, Vimentin expression was rarely detected in invasive PyMT tumor cells, further indicating that the regulatory programs inducing a trailblazer state in C3-TAg and PyMT tumors were distinct (Figure 2G, Figure 2—figure supplement 1F). The heterogeneity of Vimentin expression in C3-TAg tumors was consistent with the heterogeneity in the extent of C3-TAg organoid invasion (Figure 1A). Vimentin^pos and Vimentin^neg cells also collectively invaded together in C3-TAg primary tumors and organoids, suggesting that there was heterogeneity in collectively invading cell populations with respect to EMT state (Figure 2G, Figure 2—figure supplement 1F).

Trailblazer cells are specified by Vimentin expression in basal-like breast cancer patient tumors

We next investigated the characteristics of collectively invading cells in breast cancer patient tumors. KRT14 was expressed at a high level in basal like patient tumors compared to luminal B type breast cancer patient tumors, consistent with our results comparing the C3-TAg basal like and PyMT luminal B breast cancer GEMMs (Figure 3A). However, KRT14 expression did not correspond with the induction of collective invasion in TNBC patients, identical to the lack of association between Krt14 expression and invasion in C3-TAg primary tumor and organoids (Figure 3B, Figure 3—figure supplement 1A). In contrast, Vimentin expression was higher in the basal-like breast cancer compared to other breast cancer subtypes. Importantly, Vimentin expression was heterogeneous in TNBC tumor sections, with a trend towards more extensive Vimentin expression in regions of collective invasion (Figure 3C-D, Figure 3—figure supplement 1B). These combined results indicated that the difference in trailblazer program features in C3-TAg (basal-like) and PyMT (luminal B) tumors reflected a feature of inter-tumor heterogeneity found in human breast cancer.

Figure 3 with 1 supplement.

Metastasis is rarely detected in female mice bearing C3-TAg mammary tumors, likely because the rapid rate growth in primary tumors requires mice to be euthanized before metastases have time to develop (34). Therefore, to understand the relationship between the induction of Vimentin high cells and metastasis, we evaluated Vimentin expression in metastatic HCI-001 patient derived xenograft (PDX) TNBC tumors and their corresponding lung metastases (39). Intra- and inter-tumor heterogeneity in Vimentin expression was detected in the HCI-001 primary tumors (Figure 3E). The heterogeneity and range in Vimentin expression was similar to what was observed in C3-TAg tumors and TNBC patient tumors (Figure 2G, Figure 3D). An evaluation of the lungs revealed disseminated HCI-001 tumor cell clusters in the pulmonary vasculature (Figure 3F, top row), partially localized within the vasculature and lung tissue (Figure 3F, middle row) and fully localized within lung tissue (Figure 3F, bottom row). Notably, heterogeneous Vimentin expression was detected in each category of disseminated cell clusters, consistent with the induction of Vimentin-positive trailblazer cells metastasizing to the lungs (Figure 3F). Together, our analysis indicated that Vimentin expression specified a trailblazer state in basal-like TNBC tumors, consistent with the features of the C3-TAg trailblazer populations.

Tgfβ induces a trailblazer state in C3-TAg tumors

We next sought to define regulatory factors that activated the Vimentin specified trailblazer program. Gene set enrichment analysis (GSEA) indicated that Tgfβ signaling was more active in C3-TAg tumor organoids than non-invasive spheroids (Figure 4A, Supplementary file 2). Thus, we evaluated the contribution of Tgfβ towards regulation of a trailblazer state. A pharmacological inhibitor of Tgfβ receptor 1 (Tgfbr1), A83-01, suppressed the invasion of freshly derived C3-TAg tumor organoids (Figure 4B). 1863 cells derived in the presence of A83-01 were also less invasive than 1863T cells established from the same tumor in the absence of A83-01 (Figure 4—figure supplement 1A, Supplementary file 1). Depletion of the Tgfβ regulated transcription factor Smad3 reduced the invasion of TGFβ1 stimulated 1863 cells and 1863T cells, indicating that TGFβ regulated gene expression contributed to the observed phenotypes (Figure 4—figure supplement 1B-C). These results indicated that Tgfβ signaling was necessary to sustain a trailblazer state in C3-TAg cells.

Figure 4 with 3 supplements.

To determine if Tgfβ could initiate the conversion from an opportunist to a trailblazer state, we first established organoid lines from 3 different C3-TAg tumors in the presence of A83-01 for at least 4 weeks (Figure 4—figure supplement 2A, Supplementary file 1). These conditions yielded 1339-org, 1863-org and 1788-org tumor organoid lines that had a diminished Tgfβ activity and a reduced ability to invade relative to tumor organoids tested immediately after isolation from C3-TAg tumors (Figure 4C, Figure 4—figure supplement 2A-C, Supplementary file 1). Removal of A83-01 or SB431542, a second inhibitor of Tgfbr1, promoted tumor organoid line invasion (Figure 4D, Figure 4—figure supplement 2B-C). The invasion of the tumor organoid lines was further enhanced by exogenous TGFβ1 within 48 h (Figure 4D, Figure 4—figure supplement 2B-C, Figure 4—Videos 1-3). Persistent exposure of 1339-org cells to Tgfβ signaling for up to 22 days progressively increased the rate of collective invasion (Figure 4E, Figure 4—Videos 1-4). TGFβ signaling for at least 7 days also increased the frequency that trailblazer cells detached from collectively invading cohorts and invaded as single cells, suggesting that at least some single cell invasion reflected a loss of cohesion rather than the activation of a separate invasion program (Figure 4E, Figure 4—Video 4).

To understand how the induction of a trailblazer state by Tgfβ influenced the properties of cells during primary tumor growth, we injected control 1339-org cells cultured in A83-01 (Control) or 1339-org cells stimulated with Tgfβ for 21 days (+Tgfβ) into the mammary fat pads of immune-compromised mice. Both Control and +Tgfβ 1339-org cells initiated the growth of primary tumors (Figure 4F). Interestingly, newly isolated organoids from the Control and +Tgfβ tumors were equivalently invasive (Figure 4G). This similarity was a consequence of the induction of a trailblazer state in the Control 1339-org primary tumors, which was indicated by the new organoids derived from Control 1339-org tumors being more invasive than Control 1339-org line cells prior to tumor growth (Figure 4G). This induction of a trailblazer state was potentially due to stimulation of primary tumor cells by Tgfβ in the tumor microenvironment produced my nontumor cells, such as macrophages (40). Consistent with the equivalent invasion of the new primary tumor organoids, the tumors that developed from Control and +Tgfβ organoids contained a similar fraction of Vimentin expressing cells (Figure 4F). The equivalent invasive behavior of the tumors also corresponded with a similar number of disseminating tumor cells detected in the blood (Figure 4—figure supplement 3A). Metastatic tumor cells were not detected in the lung, indicating that the extent of dissemination was not yet sufficient to initiate colonization, consistent with metastasis being a highly inefficient process (Figure 4—figure supplement 3B). Thus, Tgfβ signaling promoted the conversion of C3-TAg cells to a trailblazer state.

Tgfβ induces multiple trailblazer states through distinct phases of re-programming

To understand how Tgfβ induced a trailblazer state, we used RNA-sequencing (RNA-seq) to analyze gene expression changes 1339-org cells after 2, 7 and 22 days after Tgfβ pathway activation. Hierarchical clustering revealed 4 different groups of genes regulated by TGFβ (Figure 5A, Supplementary file 3). If cells were undergoing binary conversion between states, the population averaged mRNA expression obtained by RNA-seq would have shown 2 groups of genes (suppressed and induced) undergoing a gradual change in expression as percentage of cells in the 2 different states shifted (41). The 4 different gene groups, each with a unique timing in their onset of change and maximal expression difference, therefore indicated that the Tgfβ stimulated C3-TAg tumor cells progressed through different stages of re-programming. Consistent with inducing multiple different states, TGFβ promoted a shift in the distribution of EpCAM expression rather than a switch between two unique EpCAM-high and EpCAM-low states, (Figure 5—figure supplement 1A). TGFβ also induced a shift in the distribution of Vimentin expression intensity in organoids at day 7 and day 22, as opposed to a switch between two distinct Vimentin^neg and Vimentin^pos states (Figure 5—figure supplement 1B).

Figure 5 with 1 supplement.

The gene groups identified by RNA-seq were associated with distinct combinations of biological processes (Figure 5B, Supplementary file 3). Of note, genes in Groups 3 and 4 were increased in expression with different kinetics in response to Tgfβ and involved with ECM regulation (Figure 5B, Supplementary file 3). This indicated that there were different mechanisms for regulating genes associated with the same biological functions in response to TGFβ. Group 3 and 4 genes were also connected with unique functions, which indicated that the nature of reprogramming changed over time and was controlled by distinct mechanisms. Interestingly, genes commonly associated with EMT were included in Group 3 (Vcam1, Itgb3, Fn1) and Group 4 (Vimentin, N-cadherin, Periostin) (Supplementary file 3). EMT-TFs were also split between Group 3 (Zeb1), Group 4 (Zeb2, Twist), or not included in any group (Snail) (Figure 5C). E-cadherin (Group 1) and the mammary epithelial lineage markers p63, Krt8 and Krt14 (Group 2) were also decreased in expression on different time scales (Figure 5C). Consistent with these results, high Vimentin and Zeb1 expression and reduced E-cadherin expression was observed in D6 and 1863T trailblazer cells, indicating that these cell lines reflected features of trailblazer cells after persistent activation of Tgfβ signaling (Figure 5—figure supplement 1C-E). Previously defined trailblazer genes were in Group 3 (Wisp1) and Group 4 (Dab2, Fstl1, Itga11, Pdgfra, Lpar1) (Figure 5D). Notably, Group 3 and Group 4 genes were enriched in the genes that were highly expressed in C3-TAg organoids, which indicated that Tgfβ re-activated a program that was present in C3-TAg tumors (Figure 2A, 5E, Supplementary file 3).

The distinct expression profile of EMT-TFs, EMT effectors and ECM regulatory genes indicated that persistent Tgfβ pathway activation induced re-programming yielded multiple different trailblazer states in response to regulatory mechanisms that evolve over time. To test this possibility we immunostained 1339-org and C3-TAg tumor organoids with Vimentin and E-Cadherin antibodies, to specify different states within the spectrum of EMT re-programming. Vimentin^neg/E-cadherin^pos, Vimentin^pos/E-cadherin^pos and Vimentin^pos/E-cadherin^neg cells were all capable of leading collective invasion (Figure 5F-G). Notably, the presence of more Vimentin^pos/E-cadherin^pos and Vimentin^pos/E-cadherin^neg trailblazer cells with increasing time after Tgfβ pathway activation correlated with an increased extent of 1339-org collective invasion (Figure 5F). Our combined results suggested that there were multiple trailblazer states induced by Tgfβ and that they had different capacities to initiate collective invasion (Figure 5H).

The transcription factors Zeb1, Zeb2 and Fra1 confer a trailblazer state

To further understand how Tgfβ induced a trailblazer state, we next determined how siRNAs targeting 12 EMT regulatory transcription factors influenced the invasion of 1863T trailblazer cells. Smad3 siRNA transfection served as a positive control. The 1863T trailblazer cells test the requirements for trailblazer cell invasion after long-term activation of the Tgfβ pathway (Figure 6—figure supplement 1A). The siRNAs targeting Fra1, Zeb1 and Zeb2 suppressed 1863T invasion with a p< 0.05 (Figure 6A), consistent with Fra1, Zeb1 and Zeb2 being the EMT-TFs that were more highly expressed in C3-TAg tumor organoids (Figure 2F). Fra1, Zeb1 and Zeb2 siRNAs also suppressed invasion to a similar extent as Smad3 siRNAs (Figure 6A). Thus we prioritized Fra1, Zeb1 and Zeb2 for further investigation.

Figure 6 with 4 supplements.

The Zeb1 and Zeb2 siRNAs suppressed 1863T cluster invasion, demonstrating that Zeb1 and Zeb2 depletion influence 1863T invasion in multiple contexts (Figure 6—figure supplement 1B). Our analysis in 1863T cells determined how Zeb1 and Zeb2 depletion influenced the invasion of trailblazer cells after long-term Tgfβ activation. To determine how Zeb1 and Zeb2 depletion influenced the induction of invasion upon the initial Tgfβ pathway stimulation, we analyzed the invasion of 1863 cells that were established and propagated in the presence of the Tgfbr1 inhibitor A83-01 before experimentation (Figure 6—figure supplement 1A). Interestingly, only Zeb1 depletion suppressed TGFβ stimulated 1863 cluster invasion, whereas Zeb2 depletion had a negligible effect (Figure 6B, Figure 6—figure supplement 1A). A miR200a-3p mimic, which suppressed Zeb1 and 2 expression, inhibited 1863 and 1863T invasion, further indicating Zeb1 and 2 contribute to the induction of a trailblazer state (Figure 6C, Figure 6—figure supplement 1C). The requirement of Zeb1, but not Zeb2, for the invasion of 1863 cells upon initial TGFβ stimulation was consistent with our results showing that Zeb1 expression rapidly increased within 7 days while increased Zeb2 expression occurred after day 7 of extrinsic TGFβ exposure (Figure 5C). Thus, our results suggested that the requirement of Zeb1 and Zeb2 may be dependent on the specific trailblazer state induced by Tgfβ signaling.

Additional investigation of Fra1 revealed that Fra1 depletion by 4 individual siRNAs and 2 different siRNA pools reduced 1863T and TGFβ stimulated 1863 invasion to a similar degree (Figure 6D, Figure 6—figure supplement 2A-D). Fra1 depletion also reduced the invasion of D6 trailblazer cells and Tgfβ stimulated 1339 cells, which were derived from a different C3-TAg tumor than the 1863 and 1863T cells (Figure 6—figure supplement 2E-G). Fra1 interacts with Jun to regulate gene expression as an AP-1 complex (42). Indeed, Jun was required for 1863T and Tgfβ stimulated 1863 cell invasion, suggesting Fra1 containing AP-1 complexes regulate the expression of genes necessary for the induction of a trailblazer state (Figure 6—figure supplement 3A-C).

Exogenous Fra1 expression promotes Zeb1 and Zeb2 expression in a mouse mammary epithelial cell line, suggesting a potential connection between the functions of Fra1, Zeb1 and Zeb2 (43). However, Fra1 depletion in 1863T cells did not alter Zeb1 or Zeb2 expression, indicating that Zeb1 and Zeb2 are not regulated by endogenous Fra1 in trailblazer cells (Figure 6E). In addition, Zeb1, Zeb2 and Vimentin expression were regulated by autocrine TGFβ signaling or exogenous TGFβ treatment whereas Fra1 expression was not influenced by autocrine or extrinsic TGFβ pathway stimulation (Figure 5C, Figure 6F-G, Supplementary file 3). These results indicated that the regulation of Fra1 expression was uncoupled from Tgfβ, Zeb1 and Zeb2. Nevertheless, Fra1 was expressed at a higher level in C3-TAg tumor organoids compared to non-invasive spheroids, consistent with increased Fra1 expression contributing to the induction of a trailblazer state (Figure 2F). In addition, FRA1 expression was higher in basal like breast cancer compared to other breast cancer subtypes (Figure 6—figure supplement 4A). FRA1 expression also correlated with VIM (Vimentin) expression in breast cancer patients, suggesting that Fra1 potentially contributes to induction of a Vimentin specified trailblazer state during tumor progression in basal-like breast cancer patients (Figure 6—figure supplement 4B). Thus, our data indicated that the induction of Fra1 and a trailblazer state in C3-TAg tumors required the activity of at least one additional signaling pathway acting in parallel with Tgfβ signaling (Figure 6H). We therefore selected Fra1 for further analysis to more fully understand the global regulatory framework that empowered a trailblazer state.

Egfr pathway activity promotes Fra1 expression and a trailblazer state

We next determined how Fra1 expression was regulated in C3-TAg tumor cells. Fra1 mRNA was more highly expressed in C3-TAg tumor organoids compared to opportunist spheroids in our RNA-seq analysis, indicating that the signaling pathway responsible for the control of Fra1 was more active in the C3-TAg organoids (Figure 2F). GSEA identified multiple active KRAS regulated expression signatures that were more active in the C3-TAg tumor organoids (Figure 7— figure supplement 1A-B, Supplementary file 2). The Kras gene can be amplified, but is not frequently mutated in C3-TAg tumors, indicating that pathway activation occurred upstream of Kras (44). Analysis of the genes in the active KRAS signature showed that the Egfr ligand Hbegf was more highly expressed in C3-TAg organoids than opportunist clones (Figure 7A, Supplementary file 2). Like FRA1, HBEGF was more highly expressed in Basal type breast tumors (Figure 7—figure supplement 1C). HBEGF expression correlated with FRA1 expression in breast cancer patient tumors as well, consistent with the possibility that Hbegf promotes Fra1 expression during tumor development (Figure 7—figure supplement 1D). Further inspection found that the Egfr ligands Areg and Epgn were also expressed at a higher level in C3-TAg organoids (Figure 7A, Supplementary file 2). These combined results prompted us to investigate how the Egfr-Ras-Mek1/2-Erk1/2 signaling pathway contributed to Fra1 expression and the induction of a trailblazer state.

Figure 7 with 4 supplements.

Erlotinib (Egfr inhibitor) and trametinib (Mek1/2 inhibitor) suppressed Fra1 mRNA and protein expression in C3-TAg cell lines and organoid lines (Figure 7B-C, Figure 7—figure supplement 2A). By comparison, Jun expression was unchanged by treatment with either inhibitor (Figure 7—figure supplement 2B). In addition, erlotinib and trametinib did not inhibit the activating phosphorylation of Smad2 or Zeb1 expression, indicating that the Egfr1-Mek1/2-Erk1/2-Fra1 did not influence cell state through autocrine Tgfβ signaling (Figure 7C, Figure 7—figure supplement 2C). A83-01 suppressed Smad2 phosphorylation, but did not influence Fra1 expression and phosphorylation, further indicating that Tgfb signaling did not regulate Fra1 (Figure 7C).

Erlotinib treatment reduced the TGFβ stimulated invasion of 1339-org cells (Figure 7D). Trametinib treatment also blocked the initial stimulation of 1339-org invasion by Tgfβ, the invasion of 1339-org trailblazer cells following persistent Tgfβ treatment and the invasion D6 trailblazer cells (Figure 7E, Figure 7—figure supplement 2D-E). Moreover, extrinsic EGF induced the invasion of 1339-org, 1863-org and 1788-org cells and this enhancement was blocked by A83-01, further indicating that the integration of Tgfbr1 and Egfr pathways promoted a trailblazer state (Figure 7F, Figure 7—figure supplement 2F-G). Trametinib inhibited the invasion of C3-Tag primary tumor organoids and phosphorylated Erk1/2 was detected in collectively invading cells in C3-TAg primary tumors (Figure 7G-H). Erk1/2 activity is pulsatile in response to EGF stimulation, and these pulses promote Fra1 expression (45). Thus, persistent Erk1/2 phosphorylation was not necessarily expected in all collectively invading cells or to promote Fra1 expression. Consistent with these results, A83-01, erlotinib and trametinib suppressed the invasion of human MCF10A-Trailblazer cells and MCFDCIS cells, which harbor features of basal-like breast cancer (25, 46). This indicated that the integration of Tgfb and Egfr pathway signaling was necessary for the induction of a trailblazer state in multiple contexts (Figure 7—figure supplement 3A-C, Figure 7—figure supplement 4A-B). Together, these results indicate that the regulation of Fra1 by Egfr and Erk1/2 promotes the induction of a trailblazer EMT (Figure 7I).

A compromise between Tgfbr1 and Egfr regulated signaling programs confers the trailblazer state

To understand the contribution of the Tgfbr1 and Egfr regulated signaling pathways towards the induction of the trailblazer state, we performed RNA-seq on organoids treated with diluent, TGFβ, TGFβ+erlotinib, EGF and EGF+A83-01 (Figure 8A, Figure 8—figure supplement 1A). Polar plots were used to cluster genes based on the nature and the degree of integration between the Tgfbr1 and Egfr signaling pathways (47). This analysis revealed multiple modes of interaction between the two pathways in response to either extrinsic TGFβ or EGF and that the type of interaction (positive or negative) was associated with different biological processes (Figure 8B-C, Figure 8—figure supplement 1B-F, Supplementary files 4 and 5). Thus, the relationship was far more complex than the simple positive signal integration previously suggested by phenotypic analysis and the expression of specific canonical EMT program genes (48, 49). In addition, the gene expression changes induced by extrinsic TGFβ and EGF were associated with unique combinations of biological processes, indicating that extrinsic TGFβ and extrinsic EGF induced distinct types of trailblazer cells (Figure 8C, Figure 8—figure supplement 1F). Nevertheless, there were conserved features of pathway integration that influenced the regulation of processes associated with invasion regardless of the external stimulus. Genes involved in controlling the organization and composition of the ECM, including previously defined trailblazer genes, were predominantly induced by Tgfbr1 and not positively impacted by Egfr signaling (Figure 8C, Figure 8—figure supplement 1F, Supplementary files 4 and 5). Notably, autocrine and extrinsic Egfr activation suppressed a number of Tgfbr1 regulated ECM reorganization genes and Hallmark EMT genes, including Postn, Bmp2, Itgb3, Orai2, Mmp11 and Erg (Figure 8D, Figure 8—figure supplement 1G, Supplementary files 4 and 5), which can contribute to invasion and metastasis (50-55). By comparison, autocrine or extrinsic EGF activity alone, or in combination with Tgfbr1 signaling, induced the expression of genes that promote cell motility, such as Axl, Sema7a, Ncam1 and Itga2 (Figure 8C and E, Figure 8—figure supplement 1F and H, Supplementary files 4 and 5) and (56-59). Thus, the induction of a set of pro-motility genes corresponded with a compromise between the Tgfbr1 and Egfr signaling programs. To determine how the Egfr dependent motility genes influenced the induction of the trailblazer state, we tested the function of the receptor tyrosine kinase Axl since it can be directly regulated by Fra1 (60). The Axl inhibitor BGB324 suppressed C3-TAg organoid invasion, demonstrating that a gene positively regulated by combined Tgfbr1 and Egfr activity was necessary for trailblazer invasion (Figure 8F-G). Together, these results suggested that the induction of the trailblazer state involved a compromise between the Tgfbr1 and Egfr regulated gene expression programs that limited the expression of genes involved in ECM reorganization and metastatic potential to promote cell motility that was essential for trailblazer cell invasion (Figure 8H).

Figure 8 with 1 supplement.

Discussion

Our findings indicate that defining how cells metastasize requires an understanding of how trailblazer cells are induced in different breast tumors (Figure 9). The Tgfβ regulated trailblazer program that we have defined in the basal-like C3-TAg mammary tumor cells induced a more aggressive form of collective invasion than was detected in luminal B PyMT tumors and functioned independently of increased Krt14 and p63 expression. This indicated that the Tgfβ program was mechanistically and phenotypically distinct from a previously described Krt14-high trailblazer program (20). In addition, PyMT tumor cells rarely expressed detectable Vimentin, which suggested that the Tgfβ regulated trailblazer program that we identified was not active in PyMT tumors. Thus, our results indicate that there are at least 2 different trailblazer programs and that their activity and requirement for invasion varies between different types of mammary tumors associated with previously defined breast cancer intrinsic subtypes (Figure 9). Genes induced as part of the C3-TAg trailblazer program are more frequently detected at a high level in basal-like breast cancer patients. Vimentin expression also correlated with collective invasion in TNBC patients, which are most often classified as basal-like. Vimentin expression was detected in basal-like TNBC PDX metastases as well. In addition, genes induced as part of the Tgfβ trailblazer program have previously been shown to correlate with reduced odds of survival in TNBC, but not patients classified as ER-positive or HER2-positive (21). These combined results indicate that the Tgfβ trailblazer program is active in basal-like TNBC and influences patient outcome. Tumor subtype specific mechanisms for inducing a trailblazer state are potentially a general property of cancer progression, as molecular subtypes are a feature of many carcinomas (61-64). Further investigation is required to determine the extent to which the Tgfβ and Krt14 trailblazer programs align with existing tumor subtypes or define new categories of inter-tumor heterogeneity. Thus, our results provide new insight into approaches used to stratify patients into outcome groups.

Figure 9.

The heterogeneity in C3-TAg tumor organoid invasion indicated that the trailblazer population itself can be heterogeneous within tumors with the same primary driving genetic factors and molecular subtype (Figure 9). Previous analysis of other tumor cell models had suggested that trailblazer cells were in a single uniform state (20, 21, 24). Notably, our analysis of 1339-org cells revealed that invasive variability can be a consequence of differences in cellular responses to a single stimulus, in this case TGFβ signaling. Therefore trailblazer cell phenotypic heterogeneity within a population did not necessarily reflect differences between EMT programs that were initiated by different activating signals (65). Gene expression changes in response to TGFβ are heterogeneous and temporally regulated in other systems as well (66). However, the influence of this gene expression heterogeneity on cell phenotype has not been clear and has often been modeled as a progressive accumulation of new functions (67). The phased changes in gene expression we defined in 1339-org cells suggested how trailblazer phenotypes evolved over time in response to TGFβ activity, with early onset gene expression changes inducing trailblazer cells that were further modulated by later stage reprogramming events (Figure 9). Thus, differences between tumor cells in the duration of TGFβ pathway activation may contribute to the heterogeneity we detect in trailblazer populations. Interestingly, long-term TGFβ induced changes in gene expression that were associated with other cell lineages. This suggests that the enhanced invasion of trailblazer cells in response to persistent TGFβ stimulation may be associated with a more general lineage plasticity and not simply a re-balancing of epithelial and mesenchymal gene expression.

Our gene expression analysis and functional studies have begun to unravel the details of the transcriptional regulation of the trailblazer state. We detected an increase in Fra1, Zeb1 and Zeb2 expression in C3-TAg trailblazer cells and human breast cancer trailblazer cell lines, with Fra1 regulated by Egfr signaling and Zeb1 and Zeb2 regulated by Tgfβ. Fra1 and Zeb1 were required for the initial Tgfβ induced induction of a trailblazer invasion as well as invasion of cells after they had converted to trailblazer state. Interestingly, Zeb2, which showed a delayed induction in response to Tgfβ, was only required for the invasion of a trailblazer cell line with the characteristics of long-term Tgfβ pathway activation. These results suggest that the transcriptional regulation of trailblazer populations may be different depending on the specific trailblazer state. We also determined how siRNAs targeting 10 additional transcription factors, including the canonical EMT TFs Snail (Snai1), Slug (Snai2) and Twist1, influenced trailblazer cell invasion. The siRNAs targeting the 10 TFs did not yield a clear suppression of invasion and we did not detect a correlation between the increased expression of the TFs and the induction of a trailblazer state. However, our results do not exclude the possibility that these TFs contribute to the Tgfβ dependent induction of a trailblazer state in some context.

Our results revealed that a new compromise mode of integration between Tgfbr1 and Egfr regulated gene expression programs. The integration of parallel signaling pathways is frequently required to confer specific cell states and is typically classified from the standpoint of being a positive or negative interaction (49, 68, 69). Indeed, prior studies have suggested that Tgfβ and Egfr-Ras regulated pathways positively integrate at the point of EMT-TF expression to initiate an EMT program (48, 49, 70, 71). While the combined activity of the Tgfbr1 and Egfr signaling pathways was necessary to induce a trailblazer state, there was an equivalent degree of positive and antagonistic interaction between the programs in C3-TAg tumor cells. This finding revealed a new type of integration between the two signaling programs in which the amplitude of a subset of the gene expression changes induced by either pathway alone are compromised. The selective activation and restriction of Tgfβ and Egf dependent gene expression suggests that compromise integration is dictated by differences in the expression or activity of co-factors that modulate chromatin accessibility or DNA binding. Tgfβ can increase Egfr expression and activate Erk1/2 signaling in a cell line dependent manner (72, 73). Similarly, Egfr-Mek1/2-Erk1/2 signaling can activate Tgfβ signaling through the phosphorylation of SMAD2/3 and increased expression of Tgfβ expression depending on cell context (43, 74, 75). However, this manner of integration is unlikely contributing to the compromise integration in the C3-TAg tumor cells given that Tgfbr1 inhibition did not alter ERK1/2 phosphorylation or Fra1 expression and Egfr and Mek1/2 inhibition failed to impact Smad2/3 phosphorylation. The direct control of the activity of one pathway by the other would also be expected to yield a more uniform positive or negative integration. Defining the mechanistic basis of selective antagonism of gene expression by the Tgfβ and Egfr pathways is an important line of future investigation.

The requirement of Tgfβ and Egfr pathway activity for trailblazer cell invasion indicated that the compromise mode of integration was a feature of the trailblazer state. Egfr pathway stimulation can also promote the motility of non-trailblazer cells, raising the possibility that Egfr signaling has distinct functions in different subpopulations of invading cells, with the precise regulatory contribution dictated by Tgfβ pathway activity (76, 77). The requirement of Tgfβ and Egfr signaling for the collective invasion of MCF10A-Trailblazer and MCFDCIS cells indicates that the integration of the pathways is essential for the induction of a trailblazer state in multiple contexts. However, the Egfr ligand Epgn restricts the invasion of luminal B PyMT organoids, potentially due to the ability of Epgn to suppress the expression of EMT-TFs and ECM reorganization genes in PyMT cells (78). This result suggests that integration of Egfr signaling is not essential for the induction of Krt14 specified trailblazer cells, like those detected in PyMT tumors (20). The differences in how Egfr signaling influences invasion may be a consequence of how tumor cells obtain the ability to resist TGFβ induced growth suppression and apoptosis (79). This is accomplished either through the clonal selection of mutants that inactivate Tgfβ signaling, which does not occur in C3-TAg tumor cells, or an intrinsic reprogramming of cellular responses to TGFβ (80). The type of reprogramming to TGFβ stimulation that is selected for during tumor progression may determine whether Tgfbr1 and Egfr signaling engage in a compromise mode of integration or if there is a predominantly positive or negative interaction between the pathways.

Our findings indicate that defining the nature of Tgfβ and Egfr pathway integration is necessary to understand tumor cell properties that influence metastasis and treatment response. In C3-TAg cells in which Tgfβ and Egfr activity combined to promote a trailblazer state, Egfr also restricted the expression of genes that could endow trailblazer cells with additional invasive and metastatic properties (53, 81-84). Similarly, Tgfβ limited the Egf induced expression of chemokines that increase invasion and the recruitment of tumor promoting immune cells (85, 86). Thus, the compromise mode of integration limited the expression of genes that could promote tumor progression whereas a strictly positive integration would have no such restrictions. The type of interaction between the Tgfβ and Egfr programs also influences the rationale for the therapeutic targeting of trailblazer cells. Inhibiting the Egfr pathway with the clinically approved drugs erlotinib and trametinib suppressed C3-TAg tumor cell invasion (87-89). However, combining either drug with newly developed anti-Tgfβ antibodies (90) may have additional benefit due to the fact that the Egfr and Tgfβ pathways function in parallel and each independently regulate large cohorts of genes. By comparison, cells in which Tgfβ activates Egfr-ERK1/2 signaling or vice versa may derive little or no benefit from combination therapy (72, 73). Thus, understanding how the mode of Tgfβ and Egfr pathway integration influences invasion, the risk of metastasis and options for treatment is an important line of future investigation.

Together, our results indicate that an improved understanding of how trailblazer cells are regulated and function in different breast cancer subtypes has the potential to improve the accuracy of patient diagnosis and reveal new treatment options.

Materials and methods

Mice

The C3-TAg [FVB-Tg(C3-1-TAg)cJeg/Jeg], PyMT [FVB/N-Tg(MMTV-PyVT)634Mul/J] and NOD.Cg-Prkdcscid Il2rgtm1Wjl/SzJ (NSG) mice were purchased from The Jackson Laboratory (C3-TAg: 013591, PyMT: 002374 and NSG: 005557) and (Supplementary file 1). Mice were housed, bred and euthanized in accordance with a protocol approved by the Institutional Animal Use and Care Committee at Georgetown University (IRB# 2017-0076) and in compliance with the NIH Guide for the Care and Use of Laboratory animals. Female mice were used for all experiments.

Patient derived xenografts

Triple negative breast cancer (TNBC) PDX fragments (Supplementary file 1) were obtained from Dr. Alana Welm at University of Utah (39). A single fragment of HCI-001 PDX was transplanted into the cleared 4^th mammary fat pad of 4 weeks old NSG female mice. Tumor measurements were taken every week to monitor tumor growth. All surgical procedures and euthanasia was performed in accordance with IACUC approved protocols at Georgetown University. Tumor pieces and lungs were fixed in formalin followed by paraffin embedding and sectioning by Histology and Tissue Shared Resource (HTSR) at Georgetown University.

Tumor organoid derivation

The largest tumors from female C3-TAg and PyMT mice were minced and dissociated for up to 120 min at 37°C in a mixture of 1 mg/ml Collagenase (Worthington, LS004188), 20 U/ml DNase (Invitrogen, 18047-019), 5% Fetal Bovine Serum (FBS) (Peak Serum, PS-FB2), ITS (Lonza, 17-838Z), non-essential amino acids (Sigma, M7145) and 10 ng/ml FGF-2 (Peprotech, 100-18C) in DMEM/F-12 (Corning, 10-092-CV). Dissociated tumors were pelleted at 80 x g for 2 min and the supernatant was discarded. Tumor organoids were then rinsed up to 5 times in 5% FBS in DMEM/F-12 followed by filtering through a 250 µm tissue strainer. Organoids were then tested in invasion assays or to establish organoid lines and cell lines (supplementary file 1).

Cell and organoid culture

Cell and organoid lines were generated from C3-TAg tumors (supplementary file 1) and tested for mycoplasma (Lonza, LT07-703) prior to use and the creation of frozen stocks. Cells and organoids were routinely used within 25 passages. Culture conditions for all cell and organoid lines are summarized in Supplementary file 1. Organoids were passaged at least once per week by dissociating tumor organoids into single cell suspensions with Dispase (Sigma, SCM133) and TryPLE (Gibco, 12605-010). The cell suspensions were plated at a density of 200,000-500,000 cells per well in a 24-well ultra-low adhesion plate (Sigma, CLS3474) to form multicellular aggregates for at least 16 h. The aggregates were then embedded in Matrigel (Corning, 354230) or Cultrex (Biotechne, 3533-005-02) and overlaid with Organoid Media (Supplementary file 1) in 6-well ultra-low adhesion plates (Sigma, CLS3471). H2B:GFP retrovirus (91) and H2B:mCherry and GFP lentivirus were produced and cells were infected as described previously (77).

Orthotopic Tumor Models

For orthotopic transplantation experiments, 1339-org cells were first cultured in Organoid media (Control) or Organoid media without A83-01 and supplemented with 2 ng/ml TGFꞵ1 (+TGFꞵ1). After 21 days of treatment, 500,000 control or +TGFꞵ1 1339-org cells were injected in the right 4^th fat pad of 8-12 weeks old female NSG mice. Seven mice per group from 2 independent experiments were analyzed. Mice were euthanized and tissue collection was performed 34 days post injection. Representative portions of primary tumors and complete lungs were fixed in formalin followed by paraffin embedding and sectioning by Histology and Tissue Shared Resource (HTSR) at Georgetown University. The remaining primary tumor portions were used for new organoid derivation and analysis.

Blood collection and Circulating Tumor Cell (CTC) isolation

For CTC isolation, mice bearing Control or +TGFꞵ1 1339-org orthotopic tumors were anesthetized with isoflurane and 800 µl of blood was extracted via cardiac puncture using a 25 G needle (BD, 305122). The blood was immediately transferred to EDTA coated micro tubes (Sarstedt, 4111395105) to avoid coagulation. Blood collected from all mice from the same treatment group (Control or +TGFꞵ1) were pooled together, mixed with red blood cell (RBC) lysis buffer (Thermo Fisher, 00433357) and incubated on ice for 10 min. The samples were then spun down at 4°C for 10 min at 250 g. The supernatant was discarded and the same steps were repeated once more to completely remove the RBCs. The cell pellet was then resuspended in cold 1X PBS and plated in poly-l-lysine pre-coated 8 well chamber slides (0.4 mL per well). The slide was placed in the 37°C incubator for at least 30 min. The cells were then fixed, permeabilized, blocked and stained with anti-Krt8 antibody to detect circulating C3-TAg tumor cells.

Organoid invasion

Organoids in 30 µl of ECM were plated onto a 20 µl base layer of ECM in 8 well chamber slides (Falcon, 354108) and allowed to invade for 48 h unless otherwise indicated. The ECM was a mixture of 2.4 mg/ml rat tail collagen I (Corning, CB-40236) and 2 mg/ml of reconstituted basement membrane (Matrigel or Cultrex) unless otherwise indicated. Tumor organoids were analyzed in a base media of DMEM/F-12 and ITS in Figure 1A, D, E and G, Figure 1—figure supplement 1B, Figure 4B, C and G, Figure 5G, Figure 7G and Figure 8G, and 1% FBS + 10 ng/ml FGF2 in Figure 1B and F. A base media of DMEM/F-12 and ITS was used for the analysis of 1339-org, 1863-org and 1788-org lines. Growth factors (2 ng/ml TGFβ1; 50 ng/ml EGF) and inhibitors were added to the base media and tested as indicated in the text and Supplementary file 1. Organoids from established organoid lines (Supplementary file 1) were tested after dissociation and aggregation for at least 16 h as described for organoid passaging. Organoids were fixed and stained with Hoechst and Phalloidin. Images were acquired with a Zeiss LSM800 laser scanning confocal microscope using 10x/0.45 (Zeiss, 1 420640-9900-000) or 20x/0.8 (Zeiss, 1 420650-9902-000) objectives. To determine the area of invasion and circularity, an image mask for each organoid was generated using ImageJ (NIH) based on the Hoechst (invasion area) or Phalloidin (circularity) signal. The total area of the invading cell nuclei was determined using the “Measure” function in ImageJ. Circularity was quantified using the “Analyze Particle” function.

siRNA experiments

Cells were transfected with 50 nM of siRNA using RNAiMax transfection reagent (Invitrogen, 13778-150) for 48 h before testing. Cells in all conditions designated as “Control” were transfected with a pool of siRNAs that does not target mouse genes. The catalog numbers for each siRNA are located in Supplementary file 1.

Vertical invasion assay

Vertical invasion assays were performed and quantified as described previously using a Collagen/Basement Membrane (BM) mix (2.4 mg/ml Collagen I and 2 mg/ml Matrigel or Cultrex) (25). Cells expressing H2B:GFP (10,000-18,000/well) were transfected in 96-well plates (Thermo, 165305) for 48 h. Media was then replaced with 50 µl Collagen/BM) ECM and overlaid with fresh media (DMEM-F12, ITS, 1% FBS, 10 ng/ml FGF2) for 48 h. Cells were imaged with a 10X/0.30 objective (Zeiss, 1 420340-9901-000) on a Zeiss LSM800 laser scanning confocal microscope. Twenty-one z-slices at 10 µm intervals over a total span of 200 μm were acquired for 9 tiled x, y positions per well. At least 3 wells (27 x,y positions total) were imaged per condition in each experiment. ImageJ software (NIH) was used to process images and quantify invasion using a custom macro as described (21). Invasion was normalized by dividing the total number of cells that had migrated ≥ 40 µm above the monolayer by the number of cells in the monolayer. The normalized invasion value for each condition was divided by the normalized invasion for the control wells transfected with an siRNA pool that does not target any mouse genes to provide a relative invasion value.

Cluster invasion

Cells (10,000-18,000/well) were transfected in 96-well plates (Corning, 3904). After 24 h, the transfected cells were trypsinized and clustered in 96-well Nunclon Sphera U-bottom low adhesion plates (Thermo Scientific, 174925) overnight (1000 cells/well). Homogeneous and mosaic clusters were also generated in U-bottom Nunclon Sphera low adhesion plates. Cells then aggregated into clusters for at least 16 h before embedding in ECM. Clusters were resuspended in 30 μl Collagen I/BM mix (2.4 mg/ml Collagen I and 2 mg/ml Matrigel or Cultrex), plated on 20 μl of a base layer of Matrigel/Collagen I, overlaid with media containing 1% FBS and 2 ng/ml TGFβ1 and allowed to invade for 24 h-48 h unless stated otherwise. The spheroids were fixed and stained with Hoechst and Phalloidin and imaged with a 10x/0.45 objective (Zeiss, 1 420640-9900-000) using a Zeiss LSM800 laser scanning confocal microscope. To determine the area of invasion, Hoechst staining was used to identify and mask nuclei invading into the ECM using ImageJ. The total area of the invading nuclei was then determined using the “Measure” function in ImageJ.

Live cell imaging

Imaging was performed using a Zeiss LSM800 laser scanning confocal microscope enclosed in a 37° C chamber supplemented with humidified CO₂. Images were acquired every 30 min for 14-18 h with a 10X/0.30 objective (Zeiss, 1 420340-9901-000). At least 15 different x,y coordinates with 5 to 7 z-slices over 50-100 μm span for each condition were acquired in parallel.

Immunofluorescence and immunoblotting

Experiments and analysis were performed as described (25) using antibodies detailed in Supplementary file 1. Immunoblots were imaged using an Odyssey scanner (Licor, 9120). Immunofluorescence images were acquired using 10x/0.45 (Zeiss, 1 420640-9900-000) and 20x/0.8 (Zeiss, 1 420650-9902-000) objectives. For analysis of the relationship between Krt14 expression and invasion, organoids were classified as invasive if they had at least two protrusions with each protrusion having at least two cells invading outward from the center mass. Images were exported as TIFFs and analyzed with ImageJ. The freehand selection tool was used to define ROIs containing individual organoids for analysis. The Krt14 fluorescence for each organoid was then defined as the “Mean Gray Value”, as determined using the “Measure” function. To define Krt14 expression in cells leading invasion, >10 organoids stained with Phalloidin, Hoechst, Krt14 and Krt8 were analyzed. Seven Z slices at 10 µm intervals over a total span of 100 µm were acquired, and the leading cells of each protrusion were designated as Krt14^low or Krt14^high. Vimentin and E-cadherin expression were analyzed in 5 µm histology sections. Organoids in 8-well chamber slides (Falcon, 354108) were washed 3x with distilled water. After detaching the chamber walls from the slides, the organoid/ECM gels were gently picked up using a single edge blade (Personna, 94-120-71) and laid out on a strip of parafilm (Bemis, PM-999). The 8-well chamber walls were placed around the gels and 150-200 µl of hydroxyethyl agarose processing gel (Thermo Fisher, HG 4000-012) was poured into the well. After the hydroxyethyl agarose processing gel solidified, it was transferred into histosettes (Simport, M498-3) and stored in 70% EtOH until embedding in paraffin by the Georgetown Histology and Tissue Shared Resource. The 5 µm sections were deparaffinized and immunostained as described (21) to detect Vimentin and E-cadherin expression (Table S1). For individual organoids, Vimentin expression was defined as the “Mean Gray Value” using the “Measure” function in ImageJ. Using ImageJ, Vimentin expression was defined as the “Mean Gray Value” for individual organoids using the “Measure” function. For characterization of the trailblazer cell state, the cell leading a collectively invading group of cells was classified as a “trailblazer cell” while cells detached from the tumor organoid were classified as “single cells”. Each cell was assigned E-caderhin^pos/ Vimentin^neg, E-Cadherin^pos/Vimentin^pos or E-Caderin^neg/ Vimentin^pos status based on E-Cadherin and Vimentin expression.

Immunohistochemistry

Sample processing and staining was performed as described (25). Images were acquired using 10x/0.45 (Zeiss, 1 420640-9900-000) and 20x/0.8 (Zeiss, 1 420650-9902-000) objectives. At least 10 regions of interest (ROIs) were randomly selected for each tumor. Cell Profiler (Broad Institute) was used to define image masks based on Krt8, Krt14 or Vimentin signal. To define their expression in tumor cells, masks based on the area of overlap between the Krt14 and Krt8 and the Vimentin and Krt8 signals were determined for each ROI. The area of individual channel and overlap masks created in CellProfiler were then loaded onto ImageJ and quantified using the “Measure” function.

Quantitative real-time PCR

Total RNA was isolated using the RNeasy Kit (Qiagen, 74104) or the PureLink™ RNA Mini Kit (Invitrogen, 12183025). Reverse transcription was performed using iScript cDNA Synthesis Kit (Bio-Rad, 1708891) following the manufacturer’s protocol. Real time PCR was performed with iQ SYBR Green Supermix (Bio-Rad, 1708882) using a RealPlex 2 thermocycler (Eppendorf, 2894). All primers were designed and purchased from Integrated DNA Technologies and their sequences are listed in Supplementary file 1.

Flow Cytometry

Cells were dissociated with TrypLe (Gibco, 12605-010) into single cell suspensions, and then washed 2x with FACs buffer (PBS + 2% FBS). Cells were then incubated with anti-CD326 (EpCAM)-PE/Dazzle 594 (Biolegend, 118236) at 4°C for 30 minutes. After staining, cells were washed with FACs buffer 2x, and stained with Helix NP Blue (Biolegend, 425305) in order to exclude dead cells. Cells were then subjected to flow cytometry analysis using BD LSR Fortessa. Results were analyzed using FCSExpress 7 (De Novo Software).

RNA sequencing

For the analysis of C3-TAg tumor organoids and non-invasive spheroids in Figure 2 and Figure 2–figure supplement 1, at least 1000 C3-TAg tumor organoids or 60,000-80,000 non-invasive cells were embedded in a 500 µl ECM mixture (2.4 mg/ml Collagen I and 2 mg/ml Matrigel or Cultrex) and plated on top of 500 uL base layer of ECM in 6 cm dishes (Falcon, 353004). C3-TAg tumor organoids were allowed to invade for 48 h. Non-invasive cells were allowed to grow and form spheroids for 6 days. Experiments were performed twice using frozen stocks of tumors and cell lines. For the analysis of 1339-org responses to growth factors and inhibitors in Figure 4, Figure 8 and Figure 8–figure supplement 1, 250,000 1339-org cells were embedded in a 100 µl ECM mixture (2.4 mg/ml Collagen I and 2 mg/ml Matrigel or Cultrex) and plated on top of a 200 µl ECM base layer in 6 cm dishes (Falcon, 353004). The gels were overlaid with control or treatment media. Organoid/ECM gels were collected in 15 mL conical tubes (Fisher, 339651), snap frozen and stored in −80° C freezer. Experiments were performed twice. For RNA isolation in all experiments, the gels were thawed on ice and disintegrated in 1 mL Trizol (Invitrogen, 15596026) by repeatedly passing the mixture through a 23G needle (BD, 305143). After the gels were homogenized, the Trizol-gel mixture was allowed to sit at RT for 5 min. Then 200 µL of chloroform (Sigma-Aldrich, C2432-25ML) was added. Following vigorous shaking for 15 sec, the mix was left for 2 min at RT to allow for phase separation. The mix was then centrifuged at 12,000 x g for 15 m at 4° C. The upper aqueous phase was then transferred to a new tube and RNA was then further purified using the RNeasy Kit (Qiagen, 74104) and or PureLink™ RNA Mini Kit (Invitrogen, 12183025). For Figure 2 and Figure 2–figure supplement 1, RNA sequencing, was performed by the Genomics and Epigenomics Shared Resource at Georgetown University. The cDNA libraries were generated using the Illumina TruSeq Stranded Total RNA Library Preparation Kit (Illumina, 20020597). Sequencing was performed using an Illumina NextSeq 500 at a depth of 12 million paired-end 75 base pair reads per sample. Data is available from GEO at GSE188312. Reads were quantified using Salmon with mouse reference genome (GRCm38/mm10). Counts were normalized and differential expression was calculated using DESeq2. For Figure 4, Figure 8 and Figure 8–figure supplement 1, RNA sequencing was performed by Novogene. RNA quality analysis was assessed by Bioanalyzer using the RNA Nano 6000 Assay Kit (Agilent, 5067-1511). The cDNA libraries were prepared and indexed using the NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, E7530L). Libraries were sequenced with the NovaSeq 6000 at a depth of 20 million paired-end 150 base pair reads per sample. Sequencing data available from the GEO at GSE188300 and GSE188323. Reads were mapped to the mouse reference genome (GRCm38/mm10) using HISAT2 and reads per gene determined using HTSeq. Counts were normalized and differential expression was calculated using DESeq2.

Polar plot analysis

Expression data was normalized and differentially expressed genes with a ≥2-fold change in expression with p< 0.05 after Benjamini-Hochberg (BH) correction were identified using the DESeq2 package in R (92, 93). Normalized expression data was collapsed for all genes as the average between condition replicates and were subset to only include genes that were differentially expressed between at least one of the three conditions (control, growth factor, growth factor plus inhibitor) assessed. To define the integration between the pathways, we modified a previously described approach for generating polar plots that has been used to identify dynamic trends in gene expression (47, 94). Relative mean signal changes over successive conditions, of control, growth factor, and growth factor plus inhibitor treatment (T°, T¹, and T², respectively), were transformed to an angle between −180 and +180 degrees (θ). This angle was further transformed to a scale from 0 to 360 degrees and then interpolated into a color based on the angle. Differentially expressed genes were plotted using polar coordinates (r,θ) with the radius (r), or distance from the focus, being equal to the value of largest absolute log₂ fold change between conditions and theta (θ) being equal to the angle corresponding to the trend of change and calculated using the following equations:

The intensity (I) is equivalent to the maximum expression for a gene across conditions and the ratio (R) between the minimum expression of a gene across conditions over the intensity, are used to calculate the angle of expression trend (θ). As such, genes are plotted as a point with an angle (θ) corresponding to relative changes in a gene’s expression with subsequent treatments at a distance from the focus corresponding to the largest magnitude of change between any pair of conditions in the treatment series. Three-point profile annotations of expression change (corresponding from left to right as T°, T¹, and T², respectively) are indicated on the outer edge of the polar scatter plots every 30 degrees in the color assigned to that particular θ. These profiles are accompanied by a radial histogram denoting the density of differentially changed genes across organoid conditions at particular θ of the polar scatterplot.

To summarize the meaning of genes with particular profiles in the polar scatterplots, genes were grouped by profiles at 30-degree intervals starting at 0 degrees, which denotes genes that were unchanged in expression from T° to T¹ but demonstrated a significant increase in expression following growth factor plus inhibitor treatment of the organoids (T²). Genes from each group were assessed for enrichment of Gene Ontology biological processes (GP) using the EnrichGo function of the clusterProfiler package in R (95). Processes with p< 0.01 after BH correction are included in Supplementary files 4 and 5. Values were transformed (-log₁₀) and displayed alongside polar scatterplots. Genes within each profile group of the polar scatterplot and utilized for GO enrichment were scaled and displayed as a line plot. A dotted line was added to these linear plots to denote the mean scaled expression pattern across organoid treatment conditions for all genes within a particular interval.

Breast Cancer Tissue Microarray

The construction of the TNBC tissue microarray (TMA) containing duplicate cores from 50 patients was performed by the Histopathology and Tissue Shared Resource (HTSR) at Georgetown University. Details of the methodology, patient consent, de-identification and demographics have been previously described (96). Multiplex staining of deparaffinized and rehydrated cores was performed through sequential rounds of staining, stripping and re-probing with antibodies recognizing KRT14, Vimentin and E-cadherin and OPAL fluorescent secondary antibodies (Supplementary file 1). The TMA cores were imaged with a Vectra3 Multi-Spectral Imaging microscope (Perkin Elmer) and analyzed with inForm software (Perkin Elmer). Tumor and stroma were segmented based on E-cadherin expression (tumor). Cores with folded sections or no clear tumor-stomal interface by manual inspection were excluded from analysis. Cores were classified “high invasion” if invasive strands were detected in >50% of the tumorstomal interface by manual inspection. The KRT14 and Vimentin signals were not visible when classifying invasion status.

Human Breast Cancer Gene Expression Analysis

The mRNA expression data for 996 patients in the Breast Invasive Carcinoma (The Cancer Genome Atlas (TCGA), PanCancer Atlas) was analyzed and downloaded from the cBioPortal (97, 98). Tumor subtype specific gene expression values are Z-scores relative to all samples. Correlative expression values are Log2 transformed. The mRNA expression for human breast cancer cells was derived from a prior analyses (21, 31) and available at the GEO (GSE58643 and GSE62569).

Statistical Methods

Data was analyzed with Prism 9.0.2 (Graphpad) and DESeq2, GSEA and ClusterProfiler packages in R. Data with a normal distribution, as determined by Shapiro-Wilk test, was analyzed by two tailed Student’s t-test. Data that did not pass a normality test were analyzed by Mann-Whitney U test or Kruskal Wallis test with Dunn’s Multiple comparison test. Differentially expressed genes with p< 0.05 after Benjamini-Hochberg correction were identified using the DESeq2 package in R. Fisher’s exact tests to determine the enrichment of Gene Ontology biological processes (BP) were performed using the EnrichGo function of the ClusterProfiler package in R. The specific statistical tests are indicated in the figure legends. Sample numbers are defined and indicated in the figure legends or on the figure panels. All experiments include a minimum of 2 biological replicates from separate biological preparations.

Data Availability

The RNA sequencing data generated in this study are publicly available at the GEO at GSE188300, GSE188312, GSE188323, GSE58643 and GSE62569. All other data are available within the article and its supplementary data files.

Materials Availability

Cell and organoid lines derived from C3-TAg tumors will be made available by the corresponding author upon reasonable request.

Disclosure of potential conflicts of interest

There are no potential conflicts of interest.

Authors’ contributions

A. Nasir: Conceptualization, investigation, writing-review and editing. S. Camacho: Conceptualization, investigation, writing-review and editing. A.T. McIntosh: Investigation, writing-review and editing. G. Graham: Investigation, writing-review and editing. R. Rahhal: Investigation, writing-review and editing. M.E. Huysman: Investigation. F. Alsharief: Investigation. A.T. Riegel: Resources, funding acquisition. G.W. Pearson: Conceptualization, investigation, writing-review and editing, funding acquisition.

Acknowledgements

Work was supported by NIH R01CA218670, Georgetown Women and Wine (to G.W. Pearson), NIH T32CA009686 (to S. Camacho, G. Graham) and NIH P30CA051008.

Figure 1—figure supplement 1.

Figure 1—Video 1. Live imaging of C3-TAg derived tumor organoids in organotypic culture (GFP, green, tumor cell). Fluorescent and bright field images were acquired for 14 h at 30 min intervals. The video corresponds with invasive heterogeneity in C3-TAg organoids shown in Figure 1A and Figure 1—figure supplement 1B.

Figure 1—Video 2. Live imaging of C3-TAg derived tumor dissociated organoids in organotypic culture. Bright field images were acquired for 5 h at 30 min intervals. The video corresponds with the intratumor invasive heterogeneity in C3-TAg tumor dissociated organoids shown in Figure 1—figure supplement 1B.

Figure 1—Video 3. Live imaging of homogenous clusters of C3-TAg tumor derived noninvasive clonal cells B6 (H2B:mCherry, red, nuclei). Fluorescent and bright field images were acquired for 14 h at 30 min intervals. The video corresponds with non-invasive B6 clusters shown in Figure 1B and Figure 1—figure supplement 1E.

Figure 1—Video 4. Live imaging of homogenous clusters of C3-TAg derived tumor cells. Bright field images were acquired for 14 h at 30 min intervals. The video corresponds with invasive C3-TAg tumor cells clusters shown in Figure 1B.

Figure 1—Video 5. Live imaging of heterogeneous clusters of C3-TAg derived tumor cells and non-invasive clonal cells B6 (H2B:mCherry, red, nuclei). Fluorescent and bright field images were acquired for 14 h at 30 min intervals. The video corresponds with opportunistic invasion of B6 cells induced by C3-TAg tumor cells as shown in Figure 1B.

Figure 1—Video 6. Live imaging of heterogeneous clusters of C3-TAg derived invasive (D6) and non-invasive (B6) clonal cells (H2B: GFP, green, nuclei; H2B:mCherry, red, nuclei). Fluorescent and bright field images were acquired for 14 h at 30 min intervals. The video corresponds with opportunistic invasion of B6 cells induced by D6 cells as shown in Figure 1— figure supplement 1E.

Figure 1—Video 7. Live imaging of PyMT derived tumor organoids in organotypic culture. Bright field images were acquired for 18 h at 30 min intervals. The video corresponds with less invasive PyMT organoids as shown in Figure 1F.

Figure 1—Video 8. Live imaging of C3-TAg derived tumor organoids in organotypic culture. Bright field images were acquired for 18 h at 30 min intervals. The video corresponds with more highly invasive C3-TAg organoids as shown in Figure 1F.

Figure 2—figure supplement 1.

Figure 3—figure supplement 1

Figure 4—figure supplement 1.

Figure 4—figure supplement 2.

Figure 4—figure supplement 3.

Figure 4—Video 1. Live imaging of C3-TAg derived 1339-org line organoids in serum starved media (H2B:GFP, green, nuclei). Fluorescent and bright field images were acquired for 18 h at 30 min intervals. The video corresponds with the baseline invasion of C3-TAg derived 1339-org line organoids as shown in Fig. 3C, S3C.

Figure 4—Video 2. Live imaging of C3-TAg derived 1339-org line organoids in serum starved media + 500 nM A830-01 (Tgfbr1 inhibitor) (H2B:GFP, green, nuclei). Fluorescent and bright field images were acquired for 18 h at 30 min intervals. The video corresponds with suppression of 1339-org line organoid invasion by inhibiting autocrine TGFβR1 activity as shown in Fig. S3C.

Figure 4—Video 3. Live imaging of C3-TAg derived 1339-org line organoids in presence of exogenous TGFβ1 (48 h) (H2B:GFP, green, nuclei). Fluorescent and bright field images were acquired for 18 h at 30 min intervals. The video corresponds with TGFβ1 induced invasion in 1339-org line organoids as shown in Fig. 3C, S3C.

Figure 4—Video 4. Live imaging of C3-TAg derived 1339-org line organoids in presence of persistent exogenous TGFβ1 (Day 16) (H2B:GFP, green, nuclei). Fluorescent and bright field images were acquired for 18 h at 30 min intervals. The video corresponds with progressively increased TGFβ1 induced invasion in 1339-org line organoids as shown in Fig. 3C.

Figure 5—figure supplement 1.

Figure 6—figure supplement 1.

Figure 6—figure supplement 2.

Figure 6—figure supplement 3.

Figure 6—figure supplement 4.

Figure 7—figure supplement 1.

Figure 7—figure supplement 2.

Figure 7—figure supplement 3.

Figure 7—figure supplement 4.

Figure 8—figure supplement 1.

Supplementary file 1. Description of GEM models, organoid/cell line models and reagents used in the study.

Supplementary file 2. RNA-Seq data (log2(CPM) genes x samples) from comparison between non-invasive clones and C3-TAg tumor organoids. Pathway enrichment analysis from comparison between non-invasive clones and C3-TAg tumor organoids

Supplementary file 3. Gene clusters and enriched biological processes from comparison between control and TGFβ1 treated C3-TAg 1339-org line organoids

Supplementary file 4. Gene clusters and enriched biological processes from comparison between control (T°), TGFβ1 (T¹) and TGFβ1 + Erlotinib (T²) treated C3-TAg 1339-org line organoids.

Supplementary file 5. Gene clusters and enriched biological processes from comparison between control (T°), EGF (T¹) and EGF + A8301 (T²) treated C3-TAg 1339-org line organoids.

Source Data Legends

Full annotated and raw immunoblots corresponding to Figure 1—figure supplement 1F, Figure 2—figure supplement 1B, Figure 5—figure supplement 1D, Figure 7C and Figure 7— figure supplement 2B-C.