"Longitudinal and multimodal auditing of tumor adaptation to CDK4/6 inhibitors in HR+ metastatic breast cancers”

CDK4/6 inhibitors (CDK4/6i) have transformed the treatment of hormone receptor-positive (HR+), HER2-negative (HR+) breast cancers as they are effective across all clinicopathological, age, and ethnicity subgroups for metastatic HR+ breast cancer. In metastatic ER+ breast cancer, CDK4/6i lead to strong and consistent improvement in survival across different lines of therapy. To understand how metastatic HR+ breast cancers become refractory to CDK4/6i, we have created a multimodal and longitudinal tumor atlas to investigate therapeutic adaptations in malignant cells and in the tumor immune microenvironment. This atlas is part of the NCI Cancer Moonshot Human Tumor Atlas Network and includes seven pairs of pre-and on-progression biopsies from five metastatic HR+ breast cancer patients treated with CDK4/6i. Biopsies were profiled with bulk genomics, transcriptomics, and proteomics as well as single-cell ATAC-seq and multiplex tissue imaging for spatial, single-cell resolution. These molecular datasets were then linked with detailed clinical metadata to create an atlas for understanding tumor adaptations during therapy. Analysis of our atlas datasets revealed a diverse but tractable set of tumor adaptations to CDK4/6i therapy. Malignant cells adapted to therapy via mTORC1 activation, cell cycle bypass, and increased replication stress. The tumor immune microenvironment displayed evidence of both immune activation and immune suppression, including increased PD-1 expression, features of T cell dysfunction, and CD163 + macrophage infiltration. Together, our metastatic ER+ breast cancer atlas represents a rich multimodal resource to understand tumor therapeutic adaptations to CDK4/6i therapy.


Introduction
Hormone receptor-positive (HR+) breast cancer is the predominant subtype of breast cancer and is defined by its malignant cells that have receptors for the hormones estrogen or progesterone.HR+ breast cancer accounts for 70% of all breast cancers and has a 5-year survival rate of 31% upon metastasis.Patients with metastatic HR+ breast cancer receive a multitude of different endocrine, targeted, and chemotherapies, deployed over many years (Waks and Winer 2019).Among the most common therapies used in metastatic HR+ breast cancers are cyclin dependent kinase 4/6 inhibitors (CDK4/6i), which have transformed care of this disease (Dean et al. 2010;Pernas et al. 2018).CDK4/6i are targeted molecular cytostatic therapies that inhibit cell cycle progression and proliferation in tumors by impeding activity of the CDK4/6-CCND1 axis.As front-line therapy for metastatic HR+ breast cancer, the combination of CDK4/6i and endocrine therapy (ET) nearly double progression-free survival rates compared to ET alone (Janni et al. 2017;Finn et al. 2016;Goetz et al. 2015).
Like other targeted molecular therapies, overcoming resistance to CDK4/6i is a key ongoing challenge for the cancer research community.HR+ breast cancers may be intrinsically resistant to CDK4/6i or may acquire resistance during therapy.Given the clinical significance of CDK4/6i in treating HR+ breast cancer, it is critical to understand mechanisms underlying resistance to these agents.Examination of clinical specimens via bulk DNA-sequencing (Wander et al. 2020), bulk RNA-sequencing (Zhu et al. 2022),and single-cell RNAsequencing (Griffiths et al. 2021) have identified genomic alterations and dysregulated oncogenic signaling as contributors to CDK4/6i resistance.Common aberrations include loss-of-function alterations of RB1 (Herrera-Abreu et al. 2016;Condorelli et al. 2018;O'Leary et al. 2018;Wander et al. 2020), and amplification or overexpression of CDK6 (Yang et al. 2017) and CCNE1/2 (Caldon et al. 2012;Herrera-Abreu et al. 2016;Wander et al. 2020); these aberrations mitigate efficacy of CDK4/6i by alleviating repression of cell cycle activity and proliferation.Moreover, these studies have further implicated mTOR (Xu et al. 2021;Knudsen et al. 2019;Romano et al. 2018), JNK/MAPK (Griffiths et al. 2021), and Hippo (Z.Li et al. 2018) pathway activation in CDK4/6i resistance.Collectively, upregulation of these pathways drives proliferation by promoting CDK4/6independent cell cycle progression, as well as protumoral tumor immune microenvironment (TIME) dynamics (Scirocchi et al. 2022;Zhu et al. 2022), including inflammatory programs that may also contribute to CDK4/6i resistance.
In previous work, we presented an omic and multidimensional spatial (OMS) atlas of a metastatic HR+ breast cancer patient treated with CDK4/6i and other therapies (Johnson et al. 2022).That atlas was created by linking detailed, longitudinal clinical metadata to clinical and exploratory multiscale molecular data analyses from four serial biopsies.In the present study, we substantially expanded our OMS atlas to create a large data resource that includes a cohort of 5 patients and 7 pre-and post-treatment biopsy pairs reflecting therapy with CDK4/6i and/or ET, and includes the founding patient in our original OMS atlas, as part of the NCI Cancer Moonshot Human Tumor Atlas Network (HTAN) (Mitri et al. 2018;Rozenblatt-Rosen et al. 2020).Herein, each biopsy in our atlas was profiled using a comprehensive suite of bulk and single-cell molecular assays, including some spatially resolved assays, linked with clinical metadata to inform analysis of the molecular datasets.Our atlas provides a rich resource for examining molecular changes within the tumor ecosystem, including genome, transcriptome, and proteome following exposure to CDK4/6i therapy, thus enabling investigation of therapeutic adaptations in both malignant cells and the TME.Malignant cell adaptations observed include genetic and epigenetic changes and increases in cell cycle activity.TIME changes observed include chronic inflammation and/or T cell suppression, and stromal tumor microenvironment barriers.This detailed atlas of metastatic breast cancers on CDK4/6i therapy significantly improves understanding of tumor adaptations to therapy and indicates potential therapeutic opportunities to combat resistance.All atlas raw data and derived datasets are available from the HTAN data portal at https://data.humantumoratlas.org/.

A Longitudinal and Multimodal Atlas of HR+ Metastatic Breast Cancers Treated with CDK4/6i
An atlas of metastatic HR+ breast cancers treated with a CDK4/6i was created by assaying paired biopsies with six complementary molecular assays (Figure 1a).Biopsy pairs were obtained from the same patient with one biopsy taken prior to starting CDK4/6i therapy, and a second biopsy taken after the CDK4/6i therapy had been discontinued due to disease progression.The atlas includes biopsy pairs from a cohort of five female patients who were profiled under the IRB-approved observational study "Molecular Mechanisms of Tumor Evolution and Resistance to Therapy" as part of HTAN (Mitri et al. 2018;Rozenblatt-Rosen et al. 2020).Patients received one of two CDK4/6i therapies, palbociclib or abemaciclib, and an ET.ET used with the CDK4/6i were either an aromatase inhibitor or fulvestrant.Patients remained on therapy until progression occurred based on provider determination.Duration of CDK4/6i treatment ranged from 1.5 months to 3.5 years (Figure 1b).Two patients, 9-1 and 9-14, received the CDK4/6i palbociclib (9-1P, 9-14P) as frontline therapy and then in later lines of therapy received abemaciclib monotherapy (9-1A, 9-14A) per its approved indication.Patient 9-1 received abemaciclib 1,129 days after palbociclib was stopped, and abemaciclib was taken for 348 days before resistance developed.Patient 9-14 received abemaciclib 309 days after palbociclib was stopped, and abemaciclib was taken for 238 days until the tumor became refractory.Additional pre-and on-progression biopsies were acquired for those patients, with a total of thirteen biopsies acquired from all five patients (Figure 1b).

Malignant cell adaptations to CDK4/6i therapy
In breast cancer, malignant epithelial cells take advantage of changes in their genomic, transcriptional, and proteomic activity to promote proliferation, differentiation, and survival.Evaluations of these datasets revealed heterogeneous malignant cell adaptations to CDK4/6i therapy across multiple modalities.
Because CDK4/6i is a cytostatic therapy that reduces cell cycle activity and proliferation, we investigated these aspects using a full suite of multimodal assays (Figure 1C) and observed a distinct split within the cohort.Two tumors, 9-14 and 9-15, displayed a relative decrease in G1 arrest scores, a transcriptional metric to quantify stalling at the G1-S checkpoint, after treatment, indicating an inability to suppress cell cycle progression despite CDK4/6 inhibition (Hafner et al. 2019).The other three tumors (9-1, 9-2, and 9-3) had positive G1 arrest scores, indicative of continued therapeutic cytostatic effects.In tumors with decreased G1 arrest scores, pathways related to cell cycle, mTORC1 signaling, and replication stress were significantly elevated following progression (p < 0.1, two-sided Mann-Whitney test) (Figure 2b).No substantial changes were observed in hormone signaling pathways post-therapy compared to pre-therapy indicating the CDK4/6i resistance mechanisms observed in this cohort were hormone-independent (Extended Data Figure 2).
Further investigation of the tumors from 9-14 and 9-15 revealed multiple adaptations potentially associated with increased proliferation and consequently therapeutic resistance.Gene and protein activity of CDK2 and E2F1 sharply increased between the pre-treatment and on-progression biopsies, indicating continued proliferation despite CDK4/6 inhibition.In 9-14P-the pair of biopsies taken before and after pabociclib treatment-this increase was accompanied by elevated CCND1 activity, cyclin D3 protein abundance, and p27 phosphorylation (T157), whereas 9-15 had increased CCNE1 activity (Figure 2c; Extended Data Figure 3).Taken together and consistent with other studies (Herrera-Abreu et al. 2016;Al-Qasem et al. 2022), observations identify tumors that bypass cyclin D-CDK4/6 signaling to progress through the G1-S checkpoint using non-canonical induction of cyclin D via CDK2 or hyperactivation of the cyclin E-CDK2 complex, to promote cell cycle progression.
In tumors from 9-14 and 9-15, high proliferative activity was accompanied by increased activity of PI3K/AKT/MTOR, a pathway frequently dysregulated in HR+ breast cancers and associated with CDK4/6i resistance (Yoshida et al. 2019;Knudsen et al. 2019;Occhipinti et al. 2020;Michaloglou et al. 2018;Herrera-Abreu et al. 2016).Activation of mTORC1 and its downstream substrates results in sustained cellular growth, delayed apoptosis, and continued cellular proliferation, including induction of cyclin D expression to promote cell cycle entry (Zhilin Zou et al. 2020;Ciołczyk-Wierzbicka et al. 2020).This may partially explain acquired resistance in 9-14P in which mTORC1 activation coincided with an increase in CCND1 expression and cyclin D3 protein abundance (Figure 2c).The mTORC1 pathway was especially elevated in 9-14P across all modalities, displaying increased phosphorylation of mTORC1 substrates (pS6 S235/S236; p4EBP1 S65) and enrichment of both chromatin accessibility and transcription of mTORC1 signaling (Figure 2b, d, and e).
Lastly, replication stress was observed in the highly proliferative tumors from 9-14 and 9-15, potentially as a result of genotoxic effects induced by aberrant cell cycle entry or prolonged cell cycle arrest (Fallah et al. 2021;Crozier et al. 2022).Both 9-14 and 9-15 had increased transcriptional and chromatin accessibility enrichment of the pathways REACTOME_ACTIVATION_OF_ATR_IN_RESPONSE_TO_REPLICATION_STRESS (Jassal et al. 2020) and G2M_CHECKPOINT (Liberzon et al. 2015).Activation of these pathways was accompanied by elevated Replication Protein A 32 (RPA32) phosphorylation, an indicator of replication stress (Figures 2b, c, and e).The role of replication stress as a mechanism of resistance is less established but hypothesized to allow cells with impaired G1/S checkpoint to escape cell cycle arrest by relying on the G2/M checkpoint machinery to surveil and repair DNA (Fallah et al. 2021;Scheidemann and Shajahan-Haq 2021).
In summary, multiple but not all tumors in this cohort displayed increased cell cycle activity and proliferation after progressing on CDK4/6i therapy that could not be accounted for by acquired genomic alterations.Cell cycle and proliferative activity were concordant across bulk RNA-seq, single cell and bulk proteomics, and scATAC-seq, indicating that perhaps, progression on CDK4/6i may be due to direct resistance to the cytostatic effects of CDK4/6i.For highly proliferative tumors, upregulation of cell cycle genes including CDK2, CCNE1, and CCND1 may bypass the CDK4/6 complex to alleviate G1 stalling and promote cell cycle progression.Hyperactivation of the PI3K/AKT/mTOR pathway was observed to be a potential mechanism of resistance.Replication stress was observed as a byproduct of tumor proliferation or as a mechanism to circumvent G1-S stalling.

Immune Responses on CDK4/6i Therapy
Previous studies have reported CDK4/6i stimulation of immunoregulatory pathways that promote anti-tumor immunity by recruiting cytotoxic T cells following antigen presentation, increasing cytokine and interferon signaling, and suppressing regulatory T cells (Spranger 2016;Goel et al. 2017;Scirocchi et al. 2022;Zhu et al. 2022).Longer PFS in HR+ breast cancer patients treated with combination palbociclib and letrozole has also been associated with expression of immunomodulatory genes (Zhu et al. 2022).These observations motivated us to quantify leukocyte contexture induced by CDK4/6i therapy within this metastatic HR+ breast cancer cohort.We used bulk RNA-seq and RPPA to quantify global transcriptional and proteomic changes and multiplexed immunohistochemistry (mIHC) (Banik et al. 2020) for spatial single-cell leukocyte profiling revealing both T cell activation and suppression were observed in most tumors.This section discusses findings supporting immune activation within the TIME of this cohort, and the next section discusses findings associated with immune suppression.
Multiple biopsy pairs from the cohort displayed increased activity of immune-related signaling pathways, such as cytolysis, interferon signaling, and Jak/Stat signaling, following CDK4/6i.Cytolytic activity, a transcriptional prognostic and putative immunotherapy biomarker (Rooney et al. 2015) that quantifies presence and activity of cytotoxic immune cells, was used to assess overall immune state in bulk RNA-seq.Cytolytic activity increased in 9-1P, 9-1A, 9-3, and 9-15 after therapy and decreased in 9-2 and 9-14P.The four biopsy pairs with increased cytolytic activity had significantly elevated enrichment of immune cell signatures and pathways following progression (p < 0.1, two-sided Mann-Whitney test) (Figure 3a; Extended Data Figures 2 and 6).Of note, the cancer hallmark pathways INTERFERON_ALPHA_RESPONSE, INTERFERON_GAMMA_RESPONSE, and IL6_JAK_STAT3_SIGNALING (Liberzon et al. 2015) were enriched, as were several other JAK/STAT pathways (Jassal et al. 2020;Kanehisa et al. 2016).Consistent with RNA-seq, three biopsy pairs (9-1P, 9-1A, 9-15) contained elevated activity of the JAK/STAT proteomic signature along with elevated abundance of multiple JAK/STAT proteins and phospho-proteins (Figure 3b; Extended Data Figure 7; Supplemental Figure S1).The association between INTERFERON_ALPHA_RESPONSE, INTERFERON_GAMMA_RESPONSE, and IL6_JAK_STAT3_SIGNALING was expected as interferons activate the JAK/STAT pathway (Villarino et al. 2015).The two tumors with decreased cytolytic activity displayed neither transcriptional nor proteomic immune activation after therapy.
While global immune signaling activation in RNA-seq supported an increased immune response, the immune system's ability to clear tumor is highly dependent on leukocyte composition and functional status within the TIME (Trujillo et al. 2018;Munhoz and Postow 2016).We used mIHC, a spatial single-cell proteomics profiling assay, to survey the types, functional status, and spatial organization of leukocytes in the TIME during CDK4/6i therapy.The four biopsy pairs with high cytolytic scores displayed distinct changes in leukocyte infiltration after therapy: tumors 9-3 and 9-15 contained higher overall CD45 + cell density, whereas 9-1P and 9-1A had decreased CD45 + cell density after therapy (Figure 3c).In tumors 9-3 and 9-15, there was a large increase in density of CD45 + cell density encompassing lineages frequently associated with anti-tumoral effects, including CD8 + T cells, helper T cells, dendritic cells, and B cells (Kravtsov et al. 2022;Marciscano and Anandasabapathy 2021;Sautès-Fridman et al. 2019) (Figure 3c, Supplemental Table 1).9-1A also exhibited increased density of helper T cells, B cells, and dendritic cells but these changes were minimal.Tumors from 9-3 and 9-15 contained increased helper T cell density after therapy, with helper T cells comprising over 15% of the CD45 + population in the 9-3 and 9-15 on-progression biopsies (Figure 3c, Supplemental Figure S2).Changes in dendritic cell abundance were highly heterogenous across biopsy pairs.Dendric cells play a critical role in antigen presentation and T cell recruitment (Marciscano and Anandasabapathy 2021); these increased in density for 9-3 and 9-15 (Figure 3c).Bulk RNA-seq results revealed increased enrichment of antigen presentation and chemokine pathways for the tumors with increased cytolytic scores and particularly for 9-15 (Figure 3a; Extended Data Figure 2) (Wang et al. 2019;Prabhakaran et al. 2017).In 9-2 and 9-14, there was consistency across assays where RNA-seq and RPPA revealed decreased immune activity after treatment, and mIHC showed decreased CD45 + cell infiltration and densities of immunoreactive cell populations (Figure 3b, c).
Given the important role of CD8 + T cells in anti-tumor immunity, we next investigated their cytotoxic state and spatial organization.CD8 + T cells induce apoptosis by secreting granzyme B and perforin (Raskov et al. 2021).Growing evidence has also associated the spatial organization of CD8 + T cells with therapeutic response and survival (Fu et al. 2021).In our cohort, CD8 + T cell populations displayed overall increase in cytotoxicity as evident by the increase in granzyme B (GRZB)-positive cells and cytolytic activity signature, suggestive of an anti-tumor immune response elicited by CDK4/6i (Figure 3a and e).Despite this overall increase in cytotoxicity, the TIME exhibited heterogeneous changes in CD8 + T cell density and spatial proximity to malignant cells on progression.Concordant with the Activated CD8 + T cell pathway (Tamborero et al. 2018) activity from RNA-seq, 9-3 and 9-15 contained increased CD8 + cell density and relative proportion (Figure 3c; Supplemental Figure S2).To assess the spatial heterogeneity of leukocyte composition within a biopsy, we divided biopsy tissue images into tiles (0.01mm 2 ; average 118 cells/tile).The 9-15 on-progression biopsy contained abundant CD8 + T cell infiltration that was highly co-localized with malignant tumor cells (Figure 3d  and e).9-3 exhibited an increase in overall CD8 + T cell density and maintained a high degree of co-localization between CD8 + T cells and tumor cells before and after therapy (Figure 3e).Although CD8 + T cells co-localized more frequently with tumor cells in 9-1A, 9-2, and 9-14, overall CD8 + T cell density decreased on-progression, indicating reduced infiltration (or retention) overall.The 9-2 pre-treatment biopsy had a high degree of spatial exclusion of CD8 + T cells distal to malignant tumor cells, indicative of potential physical barriers arising from stroma or extracellular matrix, preventing further infiltration (Figures 3e).
In summary, while four biopsy pairs showed immune activation in bulk transcriptomics and proteomics, only tumors 9-3 and 9-15 had sustained active anti-tumor immune responses based on single-cell spatial proteomics.Single-cell spatial proteomics revealed that these two tumors had substantial increases in abundance of anti-tumor immune cell populations during therapy with a shift in cytotoxic cell spatial organization toward an anti-tumor state.Two biopsy pairs from patient 9-1 on different CDK4/6i therapies exhibited increased global immune signaling in RNA-seq and RPPA data.However, single-cell spatial proteomics revealed decreased (9-1P) or only minor increases (9-1A) in the abundance of immunoreactive cell types.9-2 and 9-14 showed a reduction in immune activation consistent across bulk and single-cell spatial assays, and the pre-treatment biopsy for 9-2 exhibited TIME organization hypothesized to reduce cytotoxic cell efficacy.

T cell Suppressive Changes Observed on CDK4/6i Therapy
Many biopsies in the cohort-including those with strong evidence of immune activation-displayed evidence of immunosuppression.Tumors evade immune surveillance through multiple mechanisms to prevent T cell infiltration and promote immunomodulatory processes that reduce their activation.In T cell-inflamed environments, chronic inflammation can develop from prolonged cytokine signaling favoring tumor immunosuppression and progression (Zaidi 2019;X. Zhang et al. 2017;Qadir et al. 2017).Further, recent studies have linked interferon and JAK/STAT signaling with CDK4/6i resistance (Kettner et al. 2019;De Angelis et al. 2021).We used RNA-seq, RPPA, and mIHC to explore multiple aspects of immune suppression, including macrophage polarization, the PD-1/PD-L1 axis, regulatory T cell recruitment, and dysfunction/exhaustion of effector T cells.
Anti-inflammatory signaling caused by prolonged interferon signaling can shift myelomonocytic cell (monotypes and macrophages) polarization to protumoral states (Poh and Ernst 2018;Boutilier and Elsawa 2021;Goswami, Bose, and Baral 2021;Zijuan Zou et al. 2023) resulting in tumor progression (Larionova et al. 2020;Noy and Pollard 2014;E. Y. Lin et al. 2006;Tu et al. 2021).High intra-tumoral abundance of myelomonocytic cells represents a poor prognostic indicator of overall survival and progression-free survival in breast cancer (Larionova et al. 2020;Tiainen et al. 2015;Ni et al. 2019;Medrek et al. 2012;Esbona et al. 2018).Using mIHC, we quantified the changes of protumoral CD163 + macrophages in six biopsy pairs from our cohort.Two of six biopsy pairs (9-3 and 9-15) exhibited a substantial increase in CD163 + macrophage density after therapy, while all other biopsy pairs contained decreased CD163 + macrophage density (Figure 3c, 4a).CD163 + macrophages resided in closer proximity to malignant cells in three biopsy pairs (9-2, 9-14, 9-15) after therapy.The 9-2 pre-treatment biopsy, which was resistant to therapy, displayed striking TIME spatial organization (Figure 4b).While CD8 + T cells were sequestered away from malignant cells in the 9-2 pre-treatment biopsy, CD163 + macrophages were concentrated both in stromal regions proximal to malignant cells and infiltrated throughout malignant-cell dense regions as highlighted by the bimodal distance distribution (Figure 4a and  4b).
The PD-1/PD-L1 checkpoint axis is a key component of lymphocyte regulation and a strong indicator of immunosuppression in cancer where PD-L1 presence is negatively associated with PFS in HR+ breast cancers on CDK4/6i therapy (Zhu et al. 2022).PD-L1 is often expressed by CD163 + macrophages closely associated with malignant cells either via direct PD-L1 expression and by modulating malignant cell expression of PD-L1 (Pu and Ji 2022).Both 9-3 and 9-15 had sharp increases in bulk RNA-seq REACTOME_PD_1_SIGNALING pathway activity with moderate increases for 9-1P and 9-1A after therapy.Two pre-treatment biopsies (9-2 and 9-14) also displayed elevated activity for the same PD-1 Signaling pathway relative to the cohort, indicating the TIME for these tumors was in a T cell suppressive state prior to therapy (Figure 3a; Extended data Figure 2) (Jassal et al. 2020).PD-L1 protein abundance in mIHC was low in our cohort at 0-3% of cells, which is expected for HR+/HER2-breast cancer (Núñez Abad et al. 2022).However, small increases in cells expressing PD-L1 post therapy were observed across malignant cells, dendritic cells, and macrophages for 9-3 and 9-15 (Figure 4c).All biopsy pairs, except 9-1A, showed a decrease in the mean minimum distance between PD-1 + and PD-L1 + cells post therapy (Figure 4d).Reduced distance between PD-1 + /PD-L1 + cells allows for direct receptor-ligand interactions required to suppress T cell functional state.
Recruitment of regulatory T cells (FoxP3 + ), which function to suppress cytotoxic and helper T cells, can be promoted both by elevated CD163 + macrophages and PD-L1 expression (H.Zhang et al. 2023;Liu et al. 2011).Both T cell inflamed tumors (9-3 and 9-15) exhibited an increase in regulatory T cell density after therapy (Figure 3c, 4e).One mechanism regulatory T cells can employ to directly inhibit helper T cells is through contact-dependent signaling.We quantified the co-localization of regulatory and helper T cells using SpatialScore, a spatial biomarker associated with immunotherapy resistance in lymphoma (Phillips et al. 2021).
In the two biopsy pairs with increasing in regulatory T cell density (9-3, 9-15), the SpatialScore significantly increased or remained unchanged in the on-progression biopsy, whereas the pre-treatment biopsies for 9-1A, 9-2, and 9-14 were elevated in SpatialScore and decreased after therapy (Figure 4f).Direct interaction of regulatory T cells can shift CD8 + T cells into a dysfunctional state (Schmidt, Oberle, and Krammer 2012;C. Li et al. 2020).We quantified CD8 + T cell differentiation states to explore the extent of CD8 + T cell dysfunction using single-cell positivity of PD-1 and Eomesodermin (EOMES), biomarkers typically upregulated in exhausted or terminally differentiated dysfunctional CD8 + T cell populations.In T cell infiltrated biopsies, the majority of the CD8 + populations were composed of T effector naïve cells (PD-1 -/EOMES -), indicative of recent recruitment to the tumor, antigen inexperience, or overall lack of T cell priming (Figure 4g).9-3 displayed a substantial increase in its population of exhausted T cells (PD-1 + /EOMES + ) after therapy (24.7%) (Figure 4g) (J.Li et al. 2018).Thus indicating regulatory T cells may be impacting T cell priming capability while also pushing effector T cells toward exhaustion (Sawant et al. 2019;C. Li et al. 2020).
Our multimodal approach reveals a complex picture of both immune activation and immune suppression.Four of the six biopsy pairs analyzed with mIHC displayed substantial immune activation after CDK4/6i treatment, with two biopsy pairs having a highly T cell inflamed TME.Immune-active tumors had increased interferon signaling, JAK/STAT activation, and cytolytic activity, often in conjunction with increases in immunoreactive cell types such as cytotoxic T cells, helper T cells, and dendritic cells.Simultaneously, there was increased immunosuppressive environment after treatment in multiple tumors, including increases in CD163 + macrophages, PD-1/PD-L1 interactions, regulatory T cell recruitment, and T cell dysfunction.The concurrent immune activation and suppression indicates that many tumors become immunoreactive during therapy but eventually exhibit immunosuppression due to chronic inflammation.This arc from immunoreactive to immunosuppressive during therapy likely contributes to therapeutic resistance as T cell-inflamed tumors under prolonged antigen exposure led to ineffective tumor clearance (Trujillo et al. 2018;Spranger 2016).

Discussion
CDK4/6i therapy has improved outcomes for HR+ breast cancer patients with metastatic disease, but therapeutic resistance is exceedingly common.Studies that highlight genomic and transcriptomic mechanisms of resistance to CDK4/6i identify mechanisms of resistance in only a subset of tumors (Wander et al. 2020;Cristofanilli et al. 2018;Finn et al. 2020;Park et al. 2023).To better understand molecular mechanisms of tumor adaptation and resistance to CDK4/6i, we created a longitudinal and multimodal atlas of metastatic HR+ breast cancer during CDK4/6i therapy.This atlas includes seven pairs of biopsies from five patients where each pair had a pre-CDK4/6i biopsy and on-progression biopsy from the same patient.Using six complementary bulk and single-cell molecular assays, we deeply characterized patient biopsies to build multidimensional atlases capturing the transition as tumors develop therapeutic resistance on CDK4/6i therapy.This atlas enabled identification and quantification of tumor changes during therapy and identify candidate malignant cell and TIME mechanisms of resistance to CDK4/6i therapy.
Analysis of our atlas data indicates a diverse but tractable set of adaptations of the tumor and the TIME to CDK4/6i therapy (Figure 5a).In addition to genomic alterations acquired during therapy, malignant cell adaptations to therapy include mTORC1 pathway activation, cell-cycle bypass via increased cyclin and CDK activity, and increased replication stress.TIME adaptations included interferon and JAK/STAT signaling paired with T cell suppression and increased presence of PD-1/PD-L1 positive cells.Figures 5b and 5c summarize malignant cell and TIME adaptations on CDK4/6i therapy in each patient.At least two adaptations were identified for each tumor, and all adaptations were observed in at least two tumors.There was broad agreement across assays for support of therapeutic adaptations.While no two tumors shared the same adaptations to therapy in this limited dataset, individual adaptations were shared by subsets of tumors.Three of seven biopsy pairs, both of 9-14's biopsy pairs and 9-15, showed marked upregulation of proliferation after therapy across many assays, indicating abrogation of the cytostatic effects of CDK4/6i therapy.Four of seven biopsy pairs showed evidence of T cell activation in the form of increased cytolytic activity, increases in overall immune cell and cytotoxic cell densities, and a decrease in regulatory T cells.However, all tumors also revealed evidence of a switch to an immunosuppressive environment, including increases in PD-1/PD-L1 activity and CD163 + macrophage infiltration.Malignant cell adaptations were not correlated with TIME adaptations.
Many of the adaptations observed during CDK4/6i therapy indicate potential therapeutic vulnerabilities that present opportunity for exploitation through combination therapy with CDK4/6i or via follow up therapy to CDK4/6i.MTOR inhibitors could be combined with CDK4/6i to combat increased mTORC1 activity (Michaloglou et al. 2018).Patient 9-1 received this combination and responded for a period of time.Inhibitors targeting other CDKs may address overactivation of the CDK2-CCNE1 complex in malignant cells.ATR or WEE1 inhibitors may be effective in tumors that show increased replication stress (Fallah et al. 2021;Ubhi and Brown 2019).Immunotherapies that block PD-1/PD-L1 checkpoint activity may dampen immune evasion by the tumor while on CDK4/6i therapy, and, indeed, prior evidence that CDK4/6i primes tumors to respond to PD-1/PD-L1 inhibitors (Deng et al. 2018).Lastly, clinical trials focused on novel targeted immunotherapies may serve to counter the increased CD163 + macrophage density by inhibition or reprogramming (Fendl et al. 2023).
Importantly, these approaches could best be deployed in patients where the specific resistant pathway(s) is engaged as an adaptive response to the CDK4/6i.
Our study characterizes tumor adaptation to CDK4/6i at an exceptional depth making it possible to connect with and extend prior literature on molecular markers of resistance to CDK4/6i therapy.Previous studies have found genomic aberrations occur on CDK4/6i therapy in RB1, PIK3CA, in other cyclins, and in genes driving activation of the PI3K pathways.Aberrations in these genes and pathways arose as well, albeit at a lower frequency than previously observed (Wander et al. 2020;O'Leary et al. 2018;Park et al. 2023;Cristofanilli et al. 2018).Results from bulk RNA-seq analysis are concordant with prior studies that found resistance to CDK4/6i often is associated with transcriptional activation of the mTORC1, E2F targets, interferon signaling, and PD-1/PD-L1 pathways (Zhu et al. 2022;Park et al. 2023).Our observation of replication stress during CDK4/6i therapy matches previously reported findings (Park et al. 2023).
Through our use of single-cell spatial proteomics assays, we have characterized tumor changes on CDK4/6i in profound detail.Cell cycle activity and proliferation in malignant cells are associated with increased density and more spatial dispersion of proliferating malignant cells.Immune activation observed at the bulk transcriptomics and proteomics level was explored in depth to identify aspects of both immune activation and immune suppression in individual tumors.Tumors that displayed high degrees of immune activation frequently also presented signs of immunosuppression.This duality is consistent with previous observations regarding the effects of CDK4/6i on the TIME (Zhu et al. 2022;Teh and Aplin 2019), and with the progression of a T cellinflamed TIME in which prolonged inflammation leads to immunosuppression and a dysfunctional anti-tumor immune response (Trujillo et al. 2018;Spranger 2016).The molecular atlas that we have created provides an exceptionally deep picture of how metastatic HR+ breast cancer tumors adapt and become resistant to CDK4/6i therapy.
Previous clinical studies on CDK4/6i resistance have often relied on pre-treatment samples to identify biomarkers of response.The use of paired patient samples in our atlas taken before and after therapy provides a novel longitudinal characterization of acquired resistance to CDK4/6i.However, it is challenging to disentangle the therapeutic effects of CDK4/6i from potential mechanisms of resistance using these two time points.Future work prioritizing biopsy collection while the patient is responding to treatment may contribute toward a better understanding of the TIME state during the course of therapy.
There are several limitations to our study and its findings.Fully characterizing all tumor adaptations to CDK4/6i therapy is not possible with this small cohort given the large number of changes that we and others have observed on therapy.Additional complications in our analysis are that biopsies came from different anatomic locales, and biopsy timing was variable.Biopsy timing relative to the initiation or cessation of treatment is important when interpreting molecular measurements because drug clearance ranges from hours to days and is dependent on treatment, dosage, and patient.After drug clearance, tumors may exhibit rebound effects and hyperactive proliferative pathways (Johnston et al. 2023).While core needle biopsies have become the gold standard sampling procedure for diagnostics as a minimally invasive procedure, core needle biopsies can suffer from tissue distortion and fragmentation, disrupting spatial organization and analysis of a limited portion of the tumor.In addition, the highly heterogenous nature of metastatic lesions means biopsies may not represent the full complexity of disease within the patient or lesion, opening the possibility for sampling bias and limiting spatial power.
Overall, this atlas demonstrates how high depth molecular information captured across multiple scales and time can be used to further our understanding of mechanisms of therapeutic resistance.The highly individualized and heterogeneous nature of tumor adaptations in both malignant cells and the TIME observed across our tumor cohort highlights the need for detailed analysis of patient samples to optimize precision care.The data described here provide a rich resource capturing the molecular changes occurring in the tumor genome, transcriptome, and proteome and microenvironment during therapy.

Competing Interests
G.B.M. is a SAB member or Consultant: for Amphista, Astex, AstraZeneca, BlueDot, Chrysallis Biotechnology, Ellipses Pharma, GSK, ImmunoMET, Infinity, Ionis, Leapfrog Bio, Lilly, Medacorp, Nanostring, Nuvectis, PDX Pharmaceuticals, Qureator, Roche, Signalchem Lifesciences, Tarveda, Turbine, Zentalis Pharmaceuticals.G.B.M. has Stock/Options/Financial relationships with: Bluedot, Catena Pharmaceuticals, ImmunoMet, Nuvectis, SignalChem, Tarveda, and Turbine.G.B.M. has Licensed Technology: HRD assay to Myriad Genetics, DSP patents with Nanostring.G.B.M. has Sponsored research with AstraZeneca.J.W.G. has licensed technologies to Abbott Diagnostics; has ownership positions in Convergent Genomics, Health Technology Innovations, Zorro Bio, and PDX Pharmaceuticals; serves as a paid consultant to New Leaf Ventures; has received research support from Thermo Fisher Scientific (formerly FEI), Zeiss, Miltenyi Biotech, Cepheid (Danaher), Quantitative Imaging, Health Technology Innovations, and Micron Technologies; and owns stock in Abbott Diagnostics, AbbVie, Alphabet, Amazon, Amgen, Apple, General Electric, Gilead, Intel, Microsoft, Nvidia, and Z.I.M.mer Biomet.L.M.C. has received reagent support from Cell Signaling Technologies, Syndax Pharmaceuticals, Inc., ZielBio, Inc., and Hibercell, Inc.; holds sponsored research agreements with Syndax Pharmaceuticals, Hibercell, Inc., Prospect Creek Foundation, Lustgarten Foundation for Pancreatic Cancer Research, Susan G. Komen Foundation, and the National Foundation for Cancer Research; is on the Advisory Board for Carisma Therapeutics, Inc., CytomX Therapeutics, Inc., Shasqi, Kineta, Inc., Hibercell, Inc., Cell Signaling Technologies, Inc., Alkermes, Inc., Raska Pharma, Inc., NextCure, Guardian Bio, AstraZeneca Partner of Choice Network (OHSU Site Leader), Genenta Sciences, Pio Therapeutics Pty Ltd., and Lustgarten Foundation for Pancreatic Cancer Research Therapeutics Working Group, Inc. A.C.A. is an author on licensed patents that cover one or more aspects of the single-cell ATAC-seq technologies utilized in this study.P.K.S. is a co-founder and member of the BOD of Glencoe Software, a member of the BOD for Applied Biomath, and a member of the SAB for RareCyte, NanoString, and Montai Health; he holds equity in Glencoe, Applied Biomath, and RareCyte.P.K.S. is a consultant for Merck, and the Sorger lab has received research funding from Novartis and Merck in the past five years.Z.I.M. has a consulting/advisory role with AstraZeneca, Gilead Sciences and Daiichi Sankyo.He has also received research funding for his institution from Seattle Genetics, Novartis, AstraZeneca, Radius Health, Daiichi Sankyo, Lilly, GlaxoSmithKline, and Olema Oncology.

Human subjects
All biospecimens and data were collected under the single-center, observational study Molecular Mechanisms of Tumor Evolution and Resistance to Therapy (IRB#16113).The study was reviewed and approved by the Oregon Health & Science University (OHSU) Institutional Review Board (IRB).Eligibility criteria were determined by the enrolling physician and limited to participants at least 12 years old with localized, advanced, and/or metastatic cancer.All participants provided written informed consent to take part in the study.This study includes five female women, who at the time of consent were aged 70 (HTA9_1), 63 (HTA9_2), 57 (HTA9_3), 58 (HTA9_14), and 55 (HTA9_15).

Clinical and exploratory workflows
All biospecimens were prospectively collected by clinical research coordinators at blood draws or biopsy procedures.The cold ischemia time for all biopsy tissue specimens was between two and nine minutes starting from removal of tissue from the patient.Biospecimens were either labeled for clinical testing per institutional guidelines or deidentified for research use.All biospecimens were tracked and managed using a custom implementation of the LabVantage laboratory information management system.Clinical analytics were performed in CLIA-certified, CAP-accredited laboratories.Clinical metadata were acquired from the patient's medical record through a combination of manual and automated data abstraction.
Biospecimen designated for research analytics underwent serial sectioning.Biospecimens preserved in formalin-fixed paraffin-embedded (FFPE) blocks were sectioned into five micron sections and distributed for H&E, multiplex immunohistochemistry, and cyclic immunofluorescence (https://dx.doi.org/10.17504/protocols.io.8epv56ppjg1b/v6).Biospecimens preserved in OCT blocks were sectioned into 5-40 micron sections and distributed for H&E, cyclic immunofluorescence, and single cell ATAC sequencing (https://dx.doi.org/10.17504/protocols.io.261ge677dl47/v4).Exploratory analytics were performed in academic laboratories or research core facilities.Assay availability across patient biospecimens varied due to insufficient tissue or failing quality control.Only biospecimens collected immediately prior to initiation of CDK4/6i therapy or after progression/therapy discontinuation were included in this manuscript.Data from additional patient biospecimens not included in this manuscript are available through the HTAN Data Portal (https://data.humantumoratlas.org/).

Targeted DNA sequencing
Targeted DNA-sequencing was performed at the CLIA-licensed/CAP-accredited OHSU Knight Diagnostic Laboratories using the GeneTrails Solid Tumor Panel assay.DNA was extracted from tumor-rich regions, macro-dissected from FFPE tissue.Library preparation was performed using custom QiaSeq chemistry (QIAGEN) with multiplexed PCR and sequencing was performed using an Illumina NextSeq 500/550.The DNA library was generated using 9,228 custom-designed primer extension assays covering 613,343 base pairs across 125 to 234 cancer-related genes.The panel is routinely sequenced at an average read depth of >2,000 for highly sensitive detection of SNVs, short in/dels, and copy number alterations.All variants included in this work were reported out clinically.

Transcriptome sequencing
RNA-sequencing was performed at the CLIA-licensed/CAP-accredited OHSU Knight Diagnostic Laboratories using the RNA Transcriptome assay.RNA was extracted from tumor-rich regions, micro-dissected from FFPE tissue.Library preparation was constructed with the TruSeq RNA Library Prep Kit and sequencing was performed on the Illumina NextSeq500/550.

Reverse phase protein arrays
Flash frozen tumor tissue collected from core needle biopsies was used with the Reverse Phase Protein Array (RPPA) assay (Tibes et al. 2006).Proteomic profiling of protein and phosphoprotein was performed at the MD Anderson Cancer Center Functional Proteomics RPPA Core.

Single-cell combinatorial indexing assay for transposase-accessible chromatin sequencing
Library preparation and sequencing: Nuclei for single-cell ATAC library preparations were performed by mincing tissue on ice followed by suspension in a Liberase reaction buffer and processed according to manufacturer's protocols.Isolated nuclei were then distributed to a 96-well plate with 1,000 nuclei per well and brought up to 8 µL with water and 4µL 2× TD buffer (Nextera XT Kit, Illumina Inc. FC-131-1024).
For sciATAC preparations 1 µL of 2.5 µM indexed tagmentation complexes prepared as in Sinnamon et al. 2019(Sinnamon et al. 2019).Tagmentation was performed at 55°C.For s3ATAC preparations 1 µL of 2.5 µM indexed s3 tagmentation complexes (ScaleBio Custom Order) designed using the specifications from the publication of the s3 technology (Mulqueen et al. 2021).Tagmentation was performed at 42°C.For both workflows tagmentation was carried out for 10 minutes on an Eppendorf ThermoMixer shaking at 300 rcf and then placed on ice.Nuclei were pooled and then sorted to 12 nuclei per well using DAPI fluorescence as previously described1 into wells of a 96-well plate containing 9 µL 1× TD buffer. 1 µL of 0.1% SDS was then added per well, the plate was sealed and spun down and then incubated at 55°C for 10 minutes and then placed on ice.
Slides were incubated with primary antibodies (Supplemental Table 2) for 1 hour at room temperature or 16-17 hours at 4 °C.Primary antibody was washed off in TBST, and either anti-rat, anti-mouse, or anti-rabbit Histofine Simple Stain MAX PO horseradish peroxidase-conjugated polymer (Nichirei Biosciences, Tokyo, Japan) was applied for 30 min at room temperature, followed by AEC chromogen (Vector Laboratories, Burlingame, California, USA).
Images were acquired using the Aperio ImageScope AT (Leica Biosystems) at ×20 magnification.Following acquisition, coverslips were gently removed in 1×TBST while agitating.Removal of AEC and HRP inactivation was accomplished by incubating the slides in 0.6% fresh H2O2 in methanol for 15 minutes.AEC removal and stripping of antibodies was accomplished by ethanol gradient incubation and heat-mediated antigen retrieval such as described above between cycles.After washing and protein blocking, samples were subjected to the next round of staining.Protocol available at: https://www.protocols.io/view/mihc-staining-ohsu-coussens-39lab-sop-3i6gkhe.
Primary antibody staining: Primary antibodies were diluted in 5% NGS and 1% BSA in 1x PBS (Supplemental Table 2) and applied overnight at 4 °C in a humidified chamber, covered with plastic coverslips (IHC World, IW-2601).Following overnight incubation, tissues were washed 3 × 10 min in 1x PBS and coverslipped.
Fluorescence Signal Quenching: Following scanning, slides were soaked in a glass Coplin jar of 1x PBS until the glass coverslip could be removed without agitation, typically 10-30 minutes.Slides were placed in 10ml of freshly prepared quenching solution (20 mM sodium hydroxide (NaOH) and 3% hydrogen peroxide (H2O2) in 1x PBS) under incandescent light, for 30 minutes.Slides were removed from the chamber with forceps and washed times for 2 min in 1x PBS.The previous steps, starting with applying primary antibodies, were repeated for each round.
Scaling/Background Cohort: Expression z-scores were computed by scaling each gene relative to a background cohort containing an additional 32 HR+ metastatic breast cancer samples along with the 13 HTAN case samples (45 total).Z-scores were calculated by scaling and mean centering the log2(TPM) of each gene.

Gene Expression Signatures
Molecular Subtype Signature: Molecular subtypes were classified using the PAM50 subtype gene signature (Parker et al. 2009).A background cohort consisting of 120 breast cancer samples including the 13 samples our cohort was mean centered and used for classifying subtypes.Spearman correlations were computed between the patient samples and the pre-defined centroids of PAM50 subtypes (Parker et al. 2009).Subtypes were assigned to each sample using the highest Spearman correlation to the subtype centroid.
G1 Arrest Scores: G1 arrest scores were computed from bulk RNA-seq using a previously defined gene signature of cell cycle arrest by CDK4/6i (Hafner et al. 2019).Briefly, the change in expression z-score for 87 genes during CDK4/6i was averaged across each pair of pre-and post-treatment samples (n = 7) and multiplied by negative one.A positive score indicates G1 arrest while a negative score indicates G1 entry.
Cytolytic Activity Scores: Using bulk RNA-seq cytolytic activity scores were computed for each sample by taking the geometric mean of GZMA and PRF1 (Rooney et al. 2015).

VIPER Regulatory Network
Regulator activity was inferred with the VIPER algorithm using a TCGA BRCA ARACNe network (Alvarez et al. 2016;Lachmann et al. 2016).Briefly, VIPER infers the activity of each regulon based on the expression of its downstream targets as defined by the ARACNe network.Gene expression data used for VIPER was mean centered and scaled using a background cohort of 45 HR+ samples.

Pathway enrichment analysis
The gene set variation analysis (GSVA) algorithm was used to infer pathway activity for the 50 MSigDB Cancer Hallmark Pathways.Pathway enrichment scores were also computed with GSVA using additional curated pathways relevant to CDK4/6i adaptation from Reactome, KEGG, and BioCarta (Hänzelmann, Castelo, and Guinney 2013;Liberzon et al. 2015;Jassal et al. 2020;Kanehisa et al. 2016).These pathways include gene sets related to cell cycle regulation, mTOR activation, and JAK/STAT signaling.GSVA was performed on the log2 transformed gene expression matrix of 45 HR+ samples.The GSVA algorithm was ran using a Gaussian kernel for cumulative density function estimation and a normalized enrichment statistic of the random walk deviations.Immune cell activity was also inferred from bulk RNA-seq with the GSVA algorithm using 16 predefined gene sets of immune cell types (Tamborero et al. 2018).Three additional immune cell signatures were also included for GSVA to compute transcriptional enrichment scores representing chemokine activity, antigen presentation, and T cell inflammation (Coppola et al. 2011;Prabhakaran et al. 2017;Wang et al. 2019;Cristescu et al. 2018).

Reverse phase protein arrays
RPPA data was first normalized for protein loading for all samples within a set followed by debatching across sets to TCGA RPPA breast cancer samples using replicate-based normalization (RBN) (Akbani et al. 2014).Prior to plotting, abundance values were z-scored by scaling and mean centering relative to a background cohort consisting of the 12 cohort samples with available RPPA data plus an additional 13 HR+ metastatic breast cancer samples.These protein z-scores were used for all downstream analysis.Protein pathway activity scores were developed using previously described predictors (Johnson et al. 2022;Labrie et al. 2019).Pathway scores were computed for each sample by taking the summation of each protein predictor in a defined pathway set, times their assigned weight (positive or negative), divided by the total number of predictors in that set (Supplemental Table 3).

Single-cell combinatorial indexing assay for transposase-accessible chromatin sequencing
Sequence reads were demultiplexed using unidex (https://github.com/adeylab/unidex)matching index reads to the PCR indexes, and then for s3ATAC, the first 8 bp of Read 2 to the Tn5 index, and removing bases 9-29 from Read 2 which includes the Tn5 Mosaic End recognition sequence.Sequence data alignment and processing: FASTQ files for each biopsy were separately aligned to the GRCh38/hg38 (GCA_000001405.15)via BWA-MEM (v0.7.17-r1198-dirty) (H.Li 2013).We used the hg38 analysis set, which contains features specifically for next-generation sequencing alignment such as masking of genomic repeat arrays and a contig to redirect reads obtained from Epstein-Barr virus.Aligned reads were filtered with the scitools (Sinnamon et al. 2019) function 'bam-rmdup', which removes poorly mapped reads (quality score < 10) and collapses PCR duplicated reads on a per-barcode basis.BAM files were then split into single BAM files for each tumor and condition (pre-treatment and on-progression) as appropriate through scitools 'bam-split', and each BAM file was indexed with SAMtools (v1.10) (Danecek et al. 2021).
ArchR object generation and quality control: All single cell ATAC-seq work post-alignment was performed with the R package ArchR (v1.0.2) (Granja et al. 2021)."Arrow" files, the fundamental data storage file for ArchR projects, were generated for each biopsy's BAM file through the ArchR function createArrowFiles().Barcodes that did not meet a minimum of 3.5 log10 unique fragments or a lenient TSS enrichment cutoff of 1.5 were removed.To eliminate doublets, we ran addDoubletScores() and filterDoublets() with default settings on each Arrow file.With the exception of the pre-treatment biopsy for 9-3, which had too few barcodes for robust doublet calling, each sample had fewer than 100 doublets identified and removed.ArchR calculates perbarcode TSS enrichment and fragment counts during Arrow file generation; after our filtering was completed, patient 9-14 (Bx1 and Bx2) and 9-15 (Bx2) datasets exceeded 6.5 median TSS enrichment and 19k median fragments.Patient 9-3 sample pair were comparatively lower in quality, achieving 2.1 & 2.8 median TSS enrichment and 19k & 9.9k median fragments, respectively.Furthermore, while all on-progression datasets retained at least 1.8k passing barcodes, pre-treatment datasets for 9-3 and 9-14 had fewer than 500 each.We therefore prioritized barcode count over additional quality filters and retained a permissive TSS enrichment cutoff (1.5) for all downstream analysis.Dimensionality reduction and clustering of scATAC data: We reduced the high dimensionality of our single cell ATAC data with the iterative Latent Semantic Indexing (LSI) implementation provided by ArchR.In the standard ArchR workflow, genome-wide ATAC signal is split into 500 bp tiles for each barcode; the most accessible of these tiles are then used for the initial LSI.Preliminary clusters identified from the first round of LSI are used to identify highly variable tiles that are used for the next LSI iteration.To perform iterative LSI, we combined all scATAC-seq datasets and ran addIterativeLSI() on the top 100,000 variable tiles.Four rounds of LSI were performed at steadily increasing resolution (0.5, 1, 1.5, and 2) to help minimize batch effects.Finally, clusters of single cells were identified from the reduced dimensions via addClusters() and visualized in a 2D embedding through addUMAP() with minDist set to 0.3.
Cell typing of single cell ATAC clusters: To identify cell types in our single cell ATAC-seq data, we first identified marker genes for each cluster by providing the ArchR gene score matrix to getMarkerFeatures(). ArchR gene scores are a metric of overall chromatin accessibility for genes based on the proximity of ATAC signal within and around the gene body (Granja et al. 2021).Manual examination of significantly enriched genes (FDR < 0.05, log2 fold change > 0.5) for each cluster revealed known cell type marker genes from which we made an initial classification.We further verified marker genes by manually inspecting ATAC-seq reads for each gene in genome browser view via plotBrowserTrack().Cluster 7 was classified as endothelial based on enriched markers VWF and PECAM1 (CD31); cluster 8 was labeled immune from enrichment of PTPRC (CD45); cluster 9 was identified as fibroblast from enriched mesenchymal markers LUM, PDGFRA/B, and collagen-encoding genes.Hepatocytes (clusters 4-5) were labeled based on enriched markers ALB, TTR, and CRP.No marker genes were found for cluster 6 due to low cell count (n = 29), but we identified it as cholangiocyte based on 1) manual inspection of ATAC signal surrounding cholangiocyte markers SOX9, CFTR, KRT7 (CK7), and KRT19 (CK19), and 2) by merit of cluster 6 barcodes originating solely from liver biopsy datasets.Epithelial/tumor clusters (clusters 1, 11, and 13-15) were labeled based on their marker gene enrichment (CDH1, EPCAM, or keratins) and by merit of their barcode composition being specific to an individual patient.Clusters 10 and 12 were additionally labeled as epithelial/tumor due to evident proximity to epithelial/tumor clusters in UMAP space and corresponding patient specificity.
Gene set module scores: We calculated module scores for gene sets through the ArchR function addModuleScores().Given a set of genes, this function uses the ArchR gene score matrix to generate a quantitative measurement of chromatin accessibility for the gene group.MSigDB (Subramanian et al. 2005;Liberzon et al. 2015) "Hallmark" gene sets used for the module score analysis were obtained using the R package msigdbr (v7.5.1) (https://igordot.github.io/msigdbr/).Modules corresponding to Reactome (Matthews et al. 2009) pathways were obtained from "C2" MSigDB curated gene sets (Liberzon et al. 2011).Only tumor cells phenotyped as "tumor" were included in statistical testing.

Multiplex immunohistochemistry
Primary image processing steps were performed using the Galaxy-MCMICRO Environment (GalaxyME) (Creason et al. 2022).For each biopsy, single-channel SVS images of the whole slide for each protein assayed were aligned to a reference Hematoxylin image using PALOM (version 2022.9.1).Registration quality was validated using Avivator.Mesmer (Deepcell version 0.12.3) was used to perform nuclear and whole-cell segmentation on the resulting registered ome.tiff image using Hematoxylin, pan-cytokeratin, CD45, and CD68 marker channels.Nuclear and whole-cell masks were merged, and lone nuclei that could not be paired with a corresponding whole-cell mask were dilated by three pixels to estimate whole-cell boundaries.MCMICRO's quantification module was used to extract single cell mean marker intensities and morphological features from the refined cell mask and registered image.The quantified cell feature table was converted to AnnData format using Scimap.The single-cell intensities were log-transformed and scaled between [0-1] using Sci-kit Learn's MinMaxScaler function.Gates were selected for each marker based on the distribution of rescaled intensity values across cells in the biopsies and validated visually in the image using the Vitessce viewer (Keller et al. 2021).Cell phenotypes were called hierarchically using the assigned gates, and according to the phenotype map (Supplemental Table 1).

Cyclic immunofluorescence
Acquired images were stitched and registered across cycles based on the DAPI features using Zen Slidescan (v2.3).Each registered scene image underwent processing using a GalaxyME pipeline.Whole-cell masks were generated using Mesmer (Deepcell version 0.12.3).DAPI was used as a nuclear marker, and a maximum intensity projection of cytokeratin (CK19, CK17, CK14, CK8, CK7, CK5), CD45, and CD68 channels was used as a membrane marker.Feature extraction and intensity scaling was performed as described for mIHC.Nonhierarchical gates were selected for all markers in the panel, and downstream analyses were performed using the resulting marker-positivity calls.Gating and validation were performed as described for mIHC.

Spatial analysis
Cell density calculation: Cell densities are reported as the number of cells of a given cell type divided by the area of the tissue in mm 2 .For both imaging assays, cell centroid coordinates were converted from pixels to micrometers using the physical resolution of the platform (cycIF = 0.325 µm/pixel, mIHC = 0.5µm/pixel).Tissue area was calculated by generating a hull (alpha parameter = 0.05) using the Python package Alphashape (version 1.3.1)(https://doi.org/10.5281/zenodo.4697576)around the centroids of all segmented cells in the biopsy.The area of the hull was converted from µm 2 to mm 2 .Ripley's K of proliferating tumor cells: Cytokeratin-positive tumor cells were classified as proliferative based on positivity for the cell-cycle markers Ki67, PCNA, and pHH3.The R package spatialTIME (version 1.3.3.3)(Creed et al. 2021) was used to calculate univariate Ripley's K metrics for proliferating tumor cells for each whole slide image within a radius of 50µm.Radii between 0-100µm were tested for this analysis, and 50µm was selected to maximize inter-sample variation in the degree of clustering.This resulted in a single value per biopsy that describes the distribution of proliferating tumor cells as being randomly distributed, clustered, or dispersed in space relative to a sample-specific permuted random point process (edge_correction = 'translation', num_permutations = 100).
Grid-based spatial heterogeneity: Fishnet grids of each biopsy, with a tile size of 0.01mm 2 , were generated using the Python package GeoPandas (version 0.13.2) (http://doi.org/10.5281/zenodo.3946761).Tiles with fewer than 20 cells total were excluded from downstream analysis.Per-tile cell density distributions were calculated to assess spatial heterogeneity across the tissues.
Minimum distance calculations: The spatial_distance function from the Python package Scimap (version 1.0.0) was used to compute minimum distances between cell types of interest.
SpatialScore calculation: The SpatialScore first reported in (Phillips et al. 2021), is defined as the ratio of the distance between a helper T cell and its nearest tumor cell to the distance between a helper T cell and its nearest regulatory T cell, reported for all helper T cells in each biopsy.The Minimum distances that comprise the SpatialScore numerator and denominator were calculated using Scimap as described above.

Statistical Testing
All Mann-Whitney tests were performed in R using the wilcox.testfunction from the stats library.Two-sided tests were performed on GSVA pathway enrichment scores from RNA-seq and were not corrected for multiple testing due to low sample size.One-sided Mann-Whitney tests were performed on per-cell open chromatin accessibility enrichment scores from sciATAC-seq to evaluate concordance of pathway enrichment from bulk RNA-seq.For each tumor sample, all cells within that sample were tested against the remaining cells from all other samples using the "greater" hypothesis method.P-values were adjusted for multiple testing using the false discovery rate (FDR) method.Samples with an FDR < 0.01 were considered to be statistically significant.Only cells phenotyped as tumor were included for testing.The number of tumor cells in each sample can be seen in Extended Data Figure 5b.Two-sided Mann-Whitney tests were used between patient tumor pairs for all single-cell cycIF and mIHC data.Pairs with a p-value of < 0.01 were considered statistically significant.

Figure Legends
Figure 1: A longitudinal and multimodal atlas of five metastatic HR+ breast cancers on CDK4/6i + ET to understand how tumors adapt to this therapy combination.a) Overview of the atlas study design and data collection.Two biopsies were collected from each tumor: (1) before treatment and (2) after treatment was stopped due to disease progression.Disease progression indicates resistance to therapy.Three bulk molecular assays and three single-cell assays were applied to each biopsy, and data from these biopsies was used to identify malignant cell and TIME adaptations to therapy.Assays applied to each biopsy include bulk genomics (tDNA), transcriptomics (RNA), and proteomics (RPPA) as well as single-cell epigenomics (sciATAC) and multiplexed tissue imaging (mIHC, cycIF) for single-cell spatial proteomics.b) Patient treatment regimen and biopsy timelines.Timelines began when the pre-treatment biopsy was taken (labeled as day 0).Horizontal bars represent the duration of treatment with CDK4/6i, an estrogen receptor inhibitor (ERi), and other relevant therapies.Several biopsy pairs were collected from the same metastatic patient because the patient was treated with different CDK4/6i therapies.In total, 7 biopsy pairs were collected from 5 patients.Only therapies given between paired pre-treatment and on-progression biopsies are shown in the figure.Additional patient treatment regimens are available through the HTAN Data Portal (https://data.humantumoratlas.org/).c) Summary of assays completed on each biopsy.Clinical metadata includes CDK4/6i and ERi used for therapy, patient response, PAM50 score, biopsy site, and HTAN patient ID.
Figure 2: Multimodal profiling reveals malignant cell adaptations to CDK4/6i.a) Oncoplot of genomic alterations from high depth targeted DNA-seq.Oncoplot includes alterations from 24 genes that have been previously associated with resistance to CDK4/6i and ET and altered in at least one patient biopsy using a targeted sequencing panel of around 200 genes.Genes (rows) are grouped by pathway and function and ordered by frequency.A comprehensive set of all genomic alterations found across patient samples can be seen in extended data figure 1. b) Transcriptional gene set variation analysis (GSVA) of malignant cell pathways during CDK4/6i therapy and progression.Heatmap describes the change in pathway enrichment during each line of CDK4/6i therapy across HTAN cases.Columns are annotated with transcriptional G1 arrest scores representing the extent to which CDK4/6i prevents transcription for cell cycle progression.Statistically significant pathways between tumors with positive (n = 3) and negative (n = 4) G1 arrest scores were selected for the plot (p < 0.1, two-sided Mann-Whitney test).Pathways (rows) are grouped by category.Biopsy pairs (columns) are hierarchically clustered to show similar patterns of activity.c) Integrated heatmap of malignant cell RNA and protein changes during CDK4/6i therapy and progression.Heatmap describes the relative change of malignant cell markers during therapy using four measurements of regulation: RNA expression (RNA), regulator activity (Viper), protein abundance (Protein), and phosphoprotein (Phospho) levels.Genes and proteins included in the heatmap were selected and grouped based on GSVA results from figure 2B to highlight specific malignant cell markers of regulation and activation across bulk assays.Each value is computed by first scaling and mean centering the data relative to a background HR+ cohort and for each pair subtracting the pre-treatment biopsy from the on-progression biopsy.d) Malignant cell proteomic pathway signaling during CDK4/6i therapy and progression.Heatmap describes the change in protein pathway scores during each line of CDK4/6i therapy across HTAN cases.Pathways included in the heatmap were selected and grouped based on GSVA results from figure 2B to show concordance across RNA and protein assays.e) Violin plots showing distributions of tumor cell chromatin accessibility enrichment scores for malignant cell pathways from sciATAC-seq.Pathways were selected based on GSVA results from figure 2B to show concordance in transcriptional and epigenetic regulation across assays.For all violin plots, the mean value is shown as a point, and the first (25%) and third quartiles (75%) are represented as vertical lines within each distribution.Red asterisks represent statistically significant samples when compared to all other samples (FDR < 0.01, onesided Mann-Whitney test).f) Quantification of spatial heterogeneity of proliferating tumor cells during CDK4/6i therapy and progression using cycIF.Tumor cells were defined as positive for at least one cytokeratin marker, and were considered proliferative if positive for Ki67, PCNA, or pHH3.Top: Degree of clustering, calculated as univariate Ripley's K, for proliferative tumor cells in each biopsy within a neighborhood radius of 50 micrometers.Increasing degree of clustering suggests higher co-localization of proliferative tumor cells relative to a sample-specific permuted random point process.Ripley's K values across the cohort were all positive, suggesting that proliferative tumor cells cluster together significantly when compared to a random point process; however, the magnitude of the degree of clustering varied widely between 4172 and 24998 across the cohort.Bottom: Violin plots showing the distribution of proliferative tumor cell density across 0.01mm2 tiles within each biopsy.Mean tile densities of proliferating tumor cells varied across the cohort, with a minimum mean density of 145 cells/mm2, and a maximum of 1646 cells/mm2.For all violin plots, the mean value is shown as a point, and the first (25%) and third quartiles (75%) are represented as vertical lines within each distribution.Red asterisks represent statistical significance between patient biopsy pairs (p < 0.01, two-sided Mann-Whitney test).changes during CDK4/6i therapy and progression.Heatmap describes the relative change of TIME markers during therapy using four measurements of regulation: RNA expression (RNA), regulator activity (VIPER), protein abundance (Protein), and phosphoprotein (Phospho) levels.Genes and proteins included in the heatmap were selected and grouped based on GSVA results from figure 3A to highlight TIME markers of immune regulation and signaling across bulk assays.Each value is computed by first scaling and mean centering the data relative to a background HR+ cohort and for each pair subtracting the pre-treatment biopsy from the on-progression biopsy.c) Cell densities reported for all CD45 + immune cells and individual immune cell type populations across pre-treatment (Bx1, white) and on-progression (Bx2, dark grey) biopsies for all biopsy pairs using mIHC.Line color connecting the pre-treatment and on-progression biopsy points represents direction of change from Bx1 to Bx2.Across all biopsies, the density of all CD45 + immune cells varied between 122 cells/mm 2 and 4270 cells/mm 2 .The individual immune lineages varied in density across the cohort, ranging between 0 cells/mm 2 and 1050 cells/mm 2 .d) Pseudocolored mIHC whole slide image and select regions of interest of 9-15 Bx2 showing spatial heterogeneity of CD8 marker (magenta).For both the whole slide image and the regions of interest, the DNA channel is colored blue.The rightmost column of the regions of interest shows the pan-cytokeratin (PANCK) marker in green.On the right, analytical methods for quantifying spatial organization and heterogeneity of cell types using mIHC, including the utilization of a fishnet grid to divide the biopsy into 0.01mm2 tiles.Each tile is colored by the density of CD8 + T cells within the tile.Secondly, minimum distances are calculated between a target cell type (e.g., tumor cells) and a query cell type (e.g., CD8 + T cells) for every query cell in the sample.e) Spatial heterogeneity and cytotoxic function of CD8 + T cells during CDK4/6i therapy and progression.Left: Violin plots showing the distribution of CD8 + T cell density across 0.01mm 2 tiles from each biopsy.Mean CD8 + T cell tile densities ranged between 293 cells/mm 2 and 1348 cells/mm 2 across the cohort.Middle: Violin plot showing the per-cell distribution of minimum distances to tumor cells for all CD8 + T cells in each biopsy.Across all biopsies in the cohort, minimum distances between CD8 + T cells and tumor cells varied between 12 µm and 76 µm.Right: Bar plots showing the percentage of CD8 + T cells positive for Granzyme B. The mean percentage of Granzyme B + CD8 + T cells across the cohort was 17.4%, with a range between 2.2% and 40.13%.For all violin plots, the mean value is shown as a point, and the first (25%) and third quartiles (75%) are represented as vertical lines within each distribution.Red asterisks represent statistical significance between patient biopsy pairs (p < 0.01, two-sided Mann-Whitney test).macrophages/monocytes density across 0.01mm 2 tiles from each biopsy.Mean tile densities of CD163 + macrophages varied between 300 cells/mm 2 and 1226 cells/mm 2 across the cohort.Middle: Violin plot showing the per-cell distribution of minimum distances to tumor cells for all CD163 + macrophages/monocytes in each biopsy.Mean minimum distances between CD163 + macrophages and tumor cells varied between 12 µm and 69 µm.For all violin plots, the mean value is shown as a point, and the first (25%) and third quartiles (75%) are represented as vertical lines within each distribution.Red asterisks represent statistical significance between patient biopsy pairs (p < 0.01, two-sided Mann-Whitney test).b) Pseudocolored mIHC whole slide image and select regions of interest of 9-2 Bx1 showing spatial heterogeneity of CD8 marker (magenta) and CD163 marker (cyan).For both the whole slide image and the regions of interest, the PANCK channel is colored green.All scale bars for both the whole slide image and regions of interest are 100 µm.c) Percent of cells that are positive for PD-L1 reported for tumor cells, dendritic cells (DC), and CD163 + macrophages/monocytes across pre-treatment (Bx1, white) and on-progression (Bx2, dark grey) biopsies for all biopsy pairs using mIHC.Line color connecting the pre-treatment and on-progression biopsy points represents direction of change from Bx1 to Bx2.Across tumor cells and antigen presenting immune cell lineages, the percentage of PDL1 + cells varied between 0-3%.d) Minimum distances from PD-1 + cells to PD-L1 + cells during CDK4/6i therapy and progression reported using mIHC.For all biopsies across the cohort, the mean minimum distance between PD-1 + and PD-L1 + cells varied between 18 µm and 962 µm.For all violin plots, the mean value is shown as a point, and the first (25%) and third quartiles (75%) are represented as vertical lines within each distribution.Red asterisks represent statistical significance between patient biopsy pairs (p < 0.01, two-sided Mann-Whitney test).e) Pseudocolored mIHC images of select regions of interest from Bx1 (top) and Bx2 (bottom) of patient 9-3 showing abundance of FOXP3 (magenta), a marker for Regulatory T cells.For both images, the DNA channel is colored blue, and PANCK is green.Scale bars are 100 micrometers.f) SpatialScore during CDK4/6i therapy and progression.The SpatialScore is the ratio of the minimum distance between a helper T cell and a tumor cell to the minimum distance between a helper T cell and a Regulatory T cell, reported for all helper T cells in each biopsy.Larger SpatialScore values indicate more frequent co-location of helper T cells and Regulatory T cells.Across all biopsies, the mean SpatialScore across all helper T cells varied between 0.16 and 4.5.For all violin plots, the mean value is shown as a point, and the first (25%) and third quartiles (75%) are represented as vertical lines within each distribution.Red asterisks represent statistical significance between patient biopsy pairs (p < 0.01, two-sided Mann-Whitney test).g) Differentiation states of CD8 + T cells based on PD-1 and EOMES markers using mIHC.CD8 + T cells were classified as positive or negative for both PD-1 and EOMES and divided into groups based on the combination of these markers.Double negative (PD-1 -/EOMES -) are representative of naïve CD8 + T cells.Double positive (PD-1 + /EOMES + ) are suggestive of CD8 + T cell dysfunction and exhaustion.PD-1 + /EOMES -and PD-1 -/EOMES + are representative of early and late effector CD8 + T cells respectively.Bubble color shows the density of each T cell subset within the biopsy, and bubble size shows the relative abundance of each subset within the CD8 + T cell population.Left: Upstream activation of mTORC1 may promote G1/S entry and cell cycle progression despite CDK4/6i.Activation of the G2/M checkpoint results in G2/M stalling and increased tolerance to replication stress, further promoting uncontrolled proliferation and tumor progression.Right: Evidence suggests a mix of increasing anti-tumor immunity and immunosuppression in resistant biopsies.Increased JAK/STAT and interferon signaling, increased antigen presentation, and increased CD8 + T cell infiltration and cytotoxicity are indicative of an active anti-tumor immune response; however, high abundance of CD163 + macrophage/monocyte populations, increased Regulatory T cell infiltration, increased immune checkpoint signaling may be indicative of elevated immunosuppression.Spatial organization may also contribute to reducing the efficacy of the immune response; whereby, cytotoxic cells are spatially stratified away from tumor cells due to dense stroma and extracellular matrix.b) Summary of patient-level malignant cell adaptations to CDK4/6i across targeted DNA-seq (DNA), bulk RNA-seq (RNA), RPPA, sciATAC-seq (ATAC), and multiplex tissue imaging (MTI).For each mechanism and patient, the strength of evidence either supporting (magenta) or refuting (cyan) the mechanism is indicated for each available assay.Patient-assay support is only indicated for mechanisms in which there is a sufficient amount of evidence across available assays for a given patient.Indicators of support are left blank for all assays within a given mechanism for tumors that displayed no evidence of that mechanism contributing to CDK4/6i adaptation.Observations for patient 9-2, a clinical non-responder, are based on their pre-treatment biopsy.Tumors (rows) are annotated by CDK4/6i therapy and clinical response.c) Summary of patient-level TIME adaptations to CDK4/6i across assays.Strength of evidence across assays is indicated in 5b.

Figure 3 :
Figure 3: Multi-omics and mIHC reveal heterogeneous changes to the tumor-immune microenvironment suggestive of mechanisms of therapeutic resistance.a) Transcriptional gene set variation analysis (GSVA) of TIME pathways during CDK4/6i therapy and progression.Heatmap describes the change in pathway enrichment during each line of CDK4/6i therapy across HTAN cases.Columns are annotated with change in cytolytic activity score during treatment.Statistically significant pathways between tumors with increasing (n = 5) and decreasing (n = 2) cytolytic activity scores were selected for the plot (p < 0.1, two-sided Mann-Whitney test).Pathways (rows) are grouped by category.Biopsy pairs (columns) are hierarchically clustered to show similar patterns of activity.b) Integrated heatmap of TIME RNA and protein

Figure 4 :
Figure 4: Immune suppression changes observed on CDK4/6i therapy.a) Spatial heterogeneity of CD163 + macrophages/monocytes during CDK4/6i therapy and progression.Left: Violin plots showing the distribution of CD163+ macrophages/monocytes density across 0.01mm 2 tiles from each biopsy.Mean tile densities of CD163 + macrophages varied between 300 cells/mm 2 and 1226 cells/mm 2 across the cohort.Middle: Violin plot showing the per-cell distribution of minimum distances to tumor cells for all CD163 + macrophages/monocytes in each biopsy.Mean minimum distances between CD163 + macrophages and tumor cells varied between 12 µm and 69 µm.For all violin plots, the mean value is shown as a point, and the first (25%) and third quartiles (75%) are represented as vertical lines within each distribution.Red asterisks represent statistical significance between patient biopsy pairs (p < 0.01, two-sided Mann-Whitney test).b) Pseudocolored mIHC whole slide image and select regions of interest of 9-2 Bx1 showing spatial heterogeneity of CD8 marker (magenta) and CD163 marker (cyan).For both the whole slide image and the regions of interest, the PANCK channel is colored green.All scale bars for both the whole slide image and regions of interest are 100 µm.c) Percent of cells that are positive for PD-L1 reported for tumor cells, dendritic cells (DC), and CD163 + macrophages/monocytes across pre-treatment (Bx1, white) and on-progression (Bx2, dark grey) biopsies for all biopsy pairs using mIHC.Line color connecting the pre-treatment and on-progression biopsy points represents direction of change from Bx1 to Bx2.Across tumor cells and antigen presenting immune cell lineages, the percentage of PDL1 + cells varied between 0-3%.d) Minimum distances from PD-1 + cells to PD-L1 + cells during CDK4/6i therapy and progression reported using mIHC.For all biopsies across the cohort, the mean minimum distance between PD-1 + and PD-L1 + cells varied between 18 µm and 962 µm.For all violin plots, the mean value is shown as a point, and the first (25%) and third quartiles (75%) are represented as vertical lines within each distribution.Red asterisks represent statistical significance between patient biopsy pairs (p < 0.01, two-sided Mann-Whitney test).

Figure 5 :
Figure5: Malignant cell and tumor immune microenvironment adaptations to CDK4/6i.a) Diagram of malignant cell and tumor immune adaptations to CDK4/6i as suggested by multimodal profiling.Left: Upstream activation of mTORC1 may promote G1/S entry and cell cycle progression despite CDK4/6i.Activation of the G2/M checkpoint results in G2/M stalling and increased tolerance to replication stress, further promoting uncontrolled proliferation and tumor progression.Right: Evidence suggests a mix of increasing anti-tumor immunity and immunosuppression in resistant biopsies.Increased JAK/STAT and interferon signaling, increased antigen presentation, and increased CD8 + T cell infiltration and cytotoxicity are indicative of an active anti-tumor immune response; however, high abundance of CD163 + macrophage/monocyte populations, increased Regulatory T cell infiltration, increased immune checkpoint signaling may be indicative of elevated immunosuppression.Spatial organization may also contribute to reducing the efficacy of the immune response; whereby, cytotoxic cells are spatially stratified away from tumor cells due to dense stroma and extracellular matrix.b) Summary of patient-level malignant cell adaptations to CDK4/6i across targeted DNA-seq (DNA), bulk RNA-seq (RNA), RPPA, sciATAC-seq (ATAC), and multiplex tissue imaging (MTI).For each mechanism and patient, the strength of evidence either supporting (magenta) or refuting (cyan) the mechanism is indicated for each available assay.Patient-assay support is only indicated for mechanisms in which there is a sufficient amount of evidence across available assays for a given patient.Indicators of support are left blank for all assays within a given mechanism for tumors that displayed no evidence of that mechanism contributing to CDK4/6i adaptation.Observations for patient 9-2, a clinical non-responder, are based on their pre-treatment biopsy.Tumors (rows) are annotated by CDK4/6i therapy and clinical response.c) Summary of patient-level TIME adaptations to CDK4/6i across assays.Strength of evidence across assays is indicated in 5b.