Published: January 10, 2019

Introduction {#sec1}
============

Circulating tumor cells (CTCs) are those cells that depart from cancerous lesions and enter the bloodstream ([@bib5]). Although extraordinarily rare compared with blood cells and forced to strive for survival in circulation, CTCs are considered to be precursors of metastasis in various cancer types, including breast cancer ([@bib2], [@bib6]). CTCs are found in the blood of cancer patients as single CTCs and CTC clusters ([@bib21], [@bib42]), with the latter featuring a higher ability to seed metastasis ([@bib1]). However, it is unknown what drives their enhanced metastatic potential and what are the vulnerabilities of clustered CTCs.

Abnormal DNA methylation patterns, including both genome-wide hypomethylation and hypermethylation, have been associated with several human cancers ([@bib15], [@bib16], [@bib20], [@bib30], [@bib36]). Generally, these cancer-associated epigenetic modifications appear to affect distinct genomic areas, with hypomethylation favoring regulatory and repetitive elements versus hypermethylation, which is more frequent in CpG islands ([@bib15]). Both modifications have the ability to alter the expression of neighboring genes and contribute to the cancer phenotype ([@bib16], [@bib30]). For regulatory elements, loss of DNA methylation at transcription factor binding sites (TFBSs) can designate active transcription factor networks or networks primed for activation at later stages, e.g., during processes such as the derivation of induced pluripotent stem cells from differentiated cells ([@bib37]) or cancer progression ([@bib19]). Although DNA methylation analysis of primary tumors is extensively investigated ([@bib19], [@bib30]), the forces that shape the DNA methylome during metastatic dissemination are largely uncharacterized.

Here, we combine microfluidic-based CTC capture from breast cancer patients and mouse models, single-cell resolution DNA methylation and RNA expression analysis, a drug screen with 2,486 FDA-approved compounds, and functional validation studies in mouse models to gain insights into the biology and vulnerabilities of CTC clusters. Our study provides a genome-wide DNA methylation landscape of single and clustered CTCs in breast cancer, highlighting fundamental differences that affect metastasis and enabling the identification of cluster-targeting compounds with immediate clinical applicability.

Results {#sec2}
=======

Identification of Differentially Methylated Regions in CTC Clusters and Single CTCs {#sec2.1}
-----------------------------------------------------------------------------------

We first sought to identify active transcription factor networks by means of accessible TFBSs in single and clustered human breast CTCs, matched within individual liquid biopsies, through a genome-wide single-cell resolution DNA methylation analysis. To this end, blood samples were drawn from 43 patients with progressive breast cancer and processed with the Parsortix device ([@bib68]), a microfluidic technology that allows a size-based, antigen-agnostic enrichment of CTCs from unprocessed blood samples, specifically adapted to achieve a capture rate of \>97.2% for single CTCs and \>99.3% for CTC clusters, and no artificial cluster formation during sample processing ([Figures S1](#figs1){ref-type="fig"}A and S1B). Upon capture, live CTCs were stained for cell surface expression of EpCAM, HER2, and EGFR, and counterstained with antibodies against CD45 to identify contaminant leukocytes ([Figure S1](#figs1){ref-type="fig"}C). Upon staining verification, we identified matched single and clustered CTCs in 19% of the analyzed samples (8/43 patients), and a total of 18 marker-positive single CTCs and 29 marker-positive CTC clusters from four patients were individually micromanipulated and deposited in lysis buffer for single-cell resolution whole-genome bisulfite sequencing ([Table S1](#mmc1){ref-type="supplementary-material"}) ([@bib17], [@bib18]). In parallel, we isolated spontaneously generated GFP-labeled single CTCs and CTC clusters from three mouse xenograft models, including two human breast CTC-derived cell lines (BR16 and BRx50) and the human breast cancer cell line MDA-MB 231 (lung metastatic variant, referred to as LM2) ([@bib49], [@bib69]). In this setting, we individually micromanipulated 71 single CTCs and 48 CTC clusters ([Table S1](#mmc1){ref-type="supplementary-material"}) and also processed them for single-cell resolution whole-genome bisulfite sequencing ([@bib17], [@bib18]). Samples with a low coverage (\< 1,000 unique CpGs) or a low bisulfite conversion efficiency (CG/CHG/CHH \< 97%)---corresponding to 10.7% of patient-derived samples and 0.8% of xenograft-derived samples---were excluded from the analysis, resulting in a total of 89 single CTCs and 71 CTC clusters from patients and xenografts. On average, we achieved 3.68% CpG coverage for single CTCs and 5.86% CpG coverage for CTC clusters, in line with recent single-cell whole-genome bisulfite sequencing studies ([@bib17], [@bib18]) ([Figures S1](#figs1){ref-type="fig"}D and S1E; [Table S2](#mmc2){ref-type="supplementary-material"}). As expected, principal component analysis (PCA) mainly segregated CTCs based on the patient of origin or the specific xenograft model ([Figures S1](#figs1){ref-type="fig"}F and S1G). Metagene plot of CpG methylation revealed comparable methylation levels between single CTCs and CTC clusters across CpG islands, gene bodies, upstream (promoters) and downstream regions, including a drop of CpG methylation around the transcriptional start site, as expected ([Figures S1](#figs1){ref-type="fig"}H and S1I). We then specifically investigated differentially methylated regions (DMRs) between single CTCs and CTC clusters, evaluating average methylation levels in overlapping 5-kb windows, as previously established for single-cell DMR analysis ([@bib17], [@bib18]). For patient-derived CTCs, with this approach we identified 3,347 DMRs with a ≥ 80% methylation difference between single CTCs and CTC clusters. Of these, 1,305 regions were hypomethylated in CTC clusters and 2,042 were hypomethylated in single CTCs ([Figure 1](#fig1){ref-type="fig"}A). We then looked at xenograft-derived CTCs, and to evaluate a comparable number of DMRs as found in patients, we assessed overlapping regions with a ≥ 70% methylation difference between single CTCs and CTC clusters. We found a total of 1,430 DMRs, of which 909 hypomethylated in CTC clusters and 521 hypomethylated in single CTCs ([Figure 1](#fig1){ref-type="fig"}B). We then analyzed DMRs from both patient- and xenograft-derived CTCs using i-*cis*Target ([@bib22]). With this analysis, among hypomethylated regions that are specific to either single CTCs or CTC clusters, we found a significant enrichment for several TFBSs, many of which overlapped between patient- and xenograft-derived CTCs ([Figures 1](#fig1){ref-type="fig"}C and 1D), thus allowing us to define specific hypomethylated TFBSs that globally characterize either single CTCs or CTC clusters in both patients and xenografts. Integrated gene ontology (GO) and pathway analysis of global CTC cluster hypomethylated TFBSs revealed a remarkable enrichment for stemness-related transcription factors that coordinately regulate proliferation and pluripotency, including OCT4, NANOG, SOX2, and SIN3A, paralleling embryonic stem cell (ESCs) biology ([Figure 1](#fig1){ref-type="fig"}E) ([@bib27], [@bib47], [@bib51], [@bib58], [@bib62]). Differently, single CTCs featured hypomethylation of other TFBSs, including those that are occupied by MEF2C, JUN, MIXL1, and SHOX2, commonly enriched in various cancers ([@bib23], [@bib26], [@bib34], [@bib55]), yet independent of a core pluripotency network ([Figure 1](#fig1){ref-type="fig"}F) ([@bib27], [@bib28]). To gain insights into more subtle changes in DNA methylation occurring specifically within promoters, gene bodies, and super enhancer regions, we carried out hypergeometric-based gene set enrichment analysis of genomic features in xenograft-derived CTCs (displaying a higher homogeneity compared to patient-derived CTCs). Consistently, this analysis revealed hypermethylation and H3K27me3 repression of Polycomb-repressive complex 2 (PRC2) target gene promoters and gene bodies (including those for SUZ12 and EED) in CTC clusters ([Figures S1](#figs1){ref-type="fig"}J--S1L), as previously alluded to in cancer specimens with stem-like and proliferative features ([@bib9], [@bib31], [@bib35], [@bib56], [@bib66]) and mirroring ESCs biology ([@bib36]).Figure S1Whole-Genome Bisulfite Sequencing of Single CTCs and CTC Clusters, Related to [Figure 1](#fig1){ref-type="fig"}(A) CTC capture efficiency from blood spiked with BR16 or BRx50 single CTCs and CTC clusters, using the Parsortix device (n = 2 per cell line with 500 single CTCs and 150 CTC clusters). Error bars represent SEM.(B) 250 single BR16-GFP^+^ and 250 single BR16-RFP^+^ cells are spiked in blood and CTCs are enriched using the Parsortix device. Captured CTCs are of single color, revealing no artificial cluster formation during processing (n = 2). Error bars represent SEM.(C) Representative pictures of single CTCs and CTC clusters from breast cancer patients, enriched with the Parsortix microfluidic device and stained for EpCAM, HER2 and EGFR (*green*). White blood cells (WBCs) are counterstained with CD45 (*red*).(D and E) Bar graph showing the percent of CpG sites that are covered in individual CTC clusters and single CTCs from patients (D) and xenografts (E).(F and G) Principal component analysis of patient-derived (F) and xenograft-derived (G) single CTCs and CTC clusters, based on all features with p ≤ 0.05.(H and I) Metaplots showing the percentage (%) of CpG methylation at CpG islands (H) and reference genes (I) in CTC clusters (*blue* line) and single CTCs (*dotted red* line). TSS: Transcription Start Site; TES: Transcription End Site.(J--L) Hypergeometric gene set enrichment analysis of promoters (J), gene bodies (K) and super-enhancers (L) displaying ≥ 20% methylation difference (p value ≤ 0.01) in xenograft-derived CTC clusters compared to single CTCs. Gene sets with adjusted p value ≤ 0.05 are shown for promoters (J) and gene bodies (K). For super-enhancers (L), the top-20 significant gene sets with adjusted p value ≤ 0.05 are shown. Gene sets related to PRC2 activity are highlighted in red.(M) Histogram showing mapped reads in patient CTCs corresponding to a methylated cytosine (C) (red) or a thymine (T) (blue; corresponding to a bisulfite-converted,unmethylated cytosine) in representative regions that include binding sites for OCT4, SOX2, NANOG and SIN3A (shaded-orange box). n = number of CpGs covered.Figure 1Whole-Genome Bisulfite Sequencing Analysis of CTCs from Breast Cancer Patients and Xenografts(A) Heatmap showing methylation variable regions with ≥ 80% methylation difference between patient-derived CTC clusters and single CTCs (false discovery rate \[FDR\] \< 0.05).(B) Heatmap showing methylation variable regions with ≥ 70% methylation difference between xenograft-derived CTC clusters and single CTCs (FDR \< 0.05).(C and D) Normalized enrichment score (NES) representing enrichment (NES ≥ 3.4) of transcription factor binding sites (TFBSs) in CTC cluster hypomethylated regions (blue) and single CTC hypomethylated regions (red) of patients (C) or xenografts (D), identified using i-*cis*Target.(E and F) Integrated gene ontology (GO) and pathway enrichment analysis of TFBSs identified using i-*cis*Target in hypomethylated regions of both patient- and xenograft-derived CTC clusters (E) or single CTCs (F). The bars represent the percentage of genes detected per GO and pathway term with p ≤ 0.05.See also [Figure S1](#figs1){ref-type="fig"} and [Tables S1](#mmc1){ref-type="supplementary-material"} and [S2](#mmc2){ref-type="supplementary-material"}.

Thus, CTC clusters are clearly distinguishable from single CTCs based on their DNA methylation status at DMRs, where they mainly feature hypomethylation of bindings sites for stemness- and proliferation-associated TFs, such as OCT4, NANOG, SOX2, and SIN3A ([Figure S1](#figs1){ref-type="fig"}M, shown for patient CTCs), accompanied by a subtler hypermethylation of PRC2 target genes in promoters and gene bodies. This indicates that phenotypic differences in circulation affect DNA methylation dynamics.

Cluster-Associated Hypomethylated Regions Correlate with Poor Prognosis in Patients with Breast Cancer {#sec2.2}
------------------------------------------------------------------------------------------------------

We then tested whether the regions that are globally hypomethylated in CTC clusters are also hypomethylated in primary breast cancer. When analyzing bisulfite-sequencing data from The Cancer Genome Atlas (TCGA), we found 673/2,214 (30.39%) overlapping probes in 789 breast cancer patients, with 198 and 197 patients displaying either low (quantile Q1) or high (quantile Q4) methylation levels, respectively, in addition to a high correlation with genome-wide methylation levels ([Figures 2](#fig2){ref-type="fig"}A, [S2](#figs2){ref-type="fig"}A, and S2B). Progression-free survival (PFS) analysis on this subset of patients showed that low methylation levels within the regions that are hypomethylated in CTC clusters significantly correlate with a poor prognosis (p \< 0.05) ([Figure 2](#fig2){ref-type="fig"}B). In contrast, no correlation with PFS is observed for those regions that are hypomethylated in single CTCs ([Figures S2](#figs2){ref-type="fig"}C--S2F). Of note, while the majority of breast cancer patients in both quantile Q1 and Q4 belong to the hormone receptor-positive subtype (ER^+^/PR^+^), the proportion of patients with the HER2-positive subtype is significantly lower in quantile Q1 compared to quantile Q4 for both CTC cluster- and single CTC-associated hypomethylated regions ([Figures S2](#figs2){ref-type="fig"}G--S2J).Figure 2CTC Cluster Hypomethylated Regions Are Associated with a Poor Prognosis in Breast Cancer Patients(A) Percentage of DNA methylation (mean beta values) of overlapping probes identified in the TCGA breast cancer patient dataset. Patients are grouped in four quantiles Q1-Q4, depending on the mean DNA methylation percentage of CTC cluster hypomethylated regions.(B) Kaplan-Meier curve showing progression-free survival of breast cancer patients with overlapping probes in quantiles Q1 versus Q4 (top). The number of patients that progressed at each time point is shown (bottom).See also [Figure S2](#figs2){ref-type="fig"}.Figure S2TCGA Progression-free Survival Analysis of Breast Cancer Patients, Related to [Figure 2](#fig2){ref-type="fig"}(A) Overlap of patient- and xenograft-derived CTC cluster hypomethylated regions with TCGA probes in breast cancer patients.(B) Correlation plot of methylation levels across CTC cluster hypomethylated regions that are overlapping with TCGA probes and genome-wide methylation levels detected in the same samples.(C) Overlap of patient- and xenograft-derived single CTC hypomethylated regions with TCGA probes in breast cancer patients.(D) Percent of DNA methylation (mean beta values) of overlapping single CTC hypomethylated probes identified in the TCGA breast cancer patient dataset.Patients are grouped in four quantiles Q1-Q4, depending on the mean DNA methylation percentage of single CTC hypomethylated regions.(E) Correlation plot of methylation levels across single CTC hypomethylated regions that are overlapping with TCGA probes and genome-wide methylation levels detected in the same samples.(F) Kaplan-Meier curve showing progression-free survival of breast cancer patients with single CTC hypomethylated overlapping probes in quantiles Q1 versus Q4 (*top panel*). The number of patients that progressed at each time point of the progression-free survival analysis is shown (*bottom panel*).(G) Pie chart showing the percent of ER^+^/PR^+^, HER2^+^ and Triple Negative (TN) breast cancer patients with CTC cluster hypomethylated overlapping probes in quantiles Q1 versus Q4.(H) Bar graphs showing the percent of ER^+^/PR^+^, HER2^+^ and Triple Negative breast cancer patients in Q1 and Q4 quantiles for CTC cluster hypomethylated overlapping probes. ns, not significant. p = by Student's t test.(I) Pie chart showing the percent of ER^+^/PR^+^, HER2^+^ and Triple Negative (TN) breast cancer patients with single CTC hypomethylated overlapping probes in quantiles Q1 versus Q4.(J) Bar graphs showing the percent of ER^+^/PR^+^, HER2^+^ and Triple Negative breast cancer patients in Q1 and Q4 quantiles for single CTC hypomethylated overlapping probes. ns, not significant. p = by Student's t test.

Together, our patient PFS analysis revealed that CTC cluster-associated hypomethylated regions are also hypomethylated at the level of the primary tumor in a subset of breast cancer patients characterized by a poor prognosis, compared to patients that display higher methylation levels in the same regions.

Proliferation-Related Genes Are Enriched in CTC Clusters {#sec2.3}
--------------------------------------------------------

Our DNA methylation analysis led us to hypothesize that CTC clusters are characterized by active TF networks that support both stemness and proliferation. To identify whether the stemness- and proliferation-related TF networks are also transcriptionally active in CTC clusters compared to single CTCs, we performed single-cell resolution RNA sequencing analysis of 48 single CTCs and 24 CTC clusters, matched within individual liquid biopsies and isolated from six breast cancer patients with progressive metastatic disease, and of 49 single CTCs and 54 CTC clusters isolated from three xenograft models ([Figures S3](#figs3){ref-type="fig"}A and S3B; [Tables S1](#mmc1){ref-type="supplementary-material"} and [S3](#mmc3){ref-type="supplementary-material"}). First, we performed a transcriptome-wide weighted gene co-expression network analysis (WGCNA) and identified 32 co-expression modules, revealing gene groups that are co-enriched in either single CTCs or CTC clusters ([Figure 3](#fig3){ref-type="fig"}A). Particularly, with this approach we identified two co-expression modules, containing total of 1,976 genes, enriched in patient-derived CTC clusters (p \< 0.02; n = 1,544 for red; n = 432 for pink expression modules), while no co-expression modules were found to be significantly enriched in single CTCs ([Figure 3](#fig3){ref-type="fig"}B; [Table S4](#mmc4){ref-type="supplementary-material"}). Gene ontology (GO) network analysis of the red and pink modules revealed an enrichment of gene groups related to cell-cell junctions, cellular proliferation, and platelet activation among others ([Figure 3](#fig3){ref-type="fig"}C). We then repeated the same analysis with xenograft-derived CTCs. Transcriptome-wide WGCNA identified 21 co-expression modules ([Figure 3](#fig3){ref-type="fig"}D), of which four significantly enriched in xenograft-derived CTC clusters and containing a total of 8,332 genes (p \< 0.003; n = 159 for magenta, n = 337 for green, n = 753 for yellow, and n = 7,083 for turquoise expression modules), and three significantly enriched in xenograft-derived single CTCs and containing a total of 294 genes (p \< 0.005; n = 202 for gray, n = 41 for dark turquoise, and n = 51 for dark red) ([Figure 3](#fig3){ref-type="fig"}E; [Table S4](#mmc4){ref-type="supplementary-material"}). GO network analysis of the xenograft-derived CTC cluster networks revealed an enrichment of genes related to cell proliferation, DNA replication, cell-cell adhesion, and metabolic processes, among others ([Figure 3](#fig3){ref-type="fig"}F). In contrast, GO network analysis of the xenograft-derived single CTCs networks pointed to different processes that included RNA splicing, ATP metabolism, and ER stress response ([Figure S3](#figs3){ref-type="fig"}C). When specifically asking which processes were commonly found enriched in CTC clusters from both patients and xenografts, we identified genes related to both cellular proliferation and cell-cell adhesion ([Figure 3](#fig3){ref-type="fig"}G). While we previously reported upregulation of cell-cell adhesion components in CTC clusters ([@bib1]), whether cells within CTC clusters are also more proliferative compared to single CTCs is poorly understood. To directly address this point and validate our RNA sequencing findings, we stained CTCs from both patients and xenografts with the proliferation marker Ki67 and found that the percentage of Ki67-positive cells is indeed greatly increased in CTC clusters compared to matched single CTCs ([Figures S3](#figs3){ref-type="fig"}D--S3G).Figure S3RNA Sequencing of Patient- and Xenograft-Derived Single CTCs and CTC Clusters, Related to [Figure 3](#fig3){ref-type="fig"}(A and B) The plots show the number of detected features by number of cells in patient-derived (A) and xenograft-derived (B) CTCs processed with RNA sequencing.(C) Gene Ontology (GO) enrichment analysis of genes in the *Dark Red*, *Dark Turquoise* and *Grey* expression modules, significantly enriched in xenograft-derived single CTCs.(D--F) Dot plots showing the percent of Ki67-positive single CTCs and Ki67-positive CTCs within CTC clusters, detected in BR16 xenograft-derived CTCs (D), LM2 xenograft-derived CTCs (E), and patient 3-derived CTCs (F). ^∗^p \< 0.05 by Student's t test. Error bars represent SEM. ID = Internal ID.(G) Representative pictures of BR16 xenograft-derived single CTCs and CTC clusters, stained with Pan Cytokeratin (PanCK) (*green*), Ki67 (*red*) and DAPI (*blue*).Figure 3RNA Sequencing Analysis of Single CTCs and CTC Clusters from Breast Cancer Patients and Xenografts(A) Weighted gene co-expression network analysis (WGCNA) in patient-derived single CTCs and CTC clusters, showing a hierarchical clustering tree of co-expression modules. Each module corresponds to a branch, which is labeled by a distinct color shown underneath.(B) WGCNA identifies 32 modules with highly correlated gene expression patterns in patient-derived single CTCs and CTC clusters. Correlations between each module and CTC clusters or single CTCs are indicated by the intensity of red or green color, respectively. p value for each module is shown in brackets.(C) GO term analysis of transcripts enriched in the red and pink expression modules of patient-derived CTC clusters.(D) WGCNA in xenograft-derived single CTCs and CTC clusters, showing a hierarchical clustering tree of co-expression modules. Each module corresponds to a branch, labeled by a distinct color shown underneath.(E) WGCNA identifies 21 modules with highly correlated gene expression patterns in xenograft-derived single CTCs and CTC clusters. Correlations between each module and CTC clusters or single CTCs are indicated by the intensity of red or green color, respectively. p value for each module is shown in brackets.(F) GO term analysis of transcripts enriched in the magenta, green, yellow, and turquoise expression modules of xenograft-derived CTC clusters.(G) Venn diagram showing the overlap of GO terms enriched in CTC clusters from patients and xenografts.See also [Figure S3](#figs3){ref-type="fig"} and [Tables S3](#mmc3){ref-type="supplementary-material"} and [S4](#mmc4){ref-type="supplementary-material"}.

In summary, our transcriptome-wide WGCNA revealed that CTC clusters---additionally to upregulating cell-cell junction components---are also characterized by a higher proliferation rate compared to single CTCs, in line with our DNA methylation results.

Stem-Cell-Related Genes Are Enriched in CTC Clusters {#sec2.4}
----------------------------------------------------

Our transcriptome-wide analysis revealed that the most significant differences between single and clustered CTCs involve the expression of cell-cell junction components and cell cycle progression. Yet our DNA methylation analysis pointed to both proliferation- and stemness-related TF networks being accessible in CTC clusters. To specifically ask whether CTC clusters would differ from single CTCs also in regard to the expression of stemness-related genes (which may be regulated with more subtle expression changes compared to cell cycle- and cell-cell junction-related genes), we focused our analysis on 335 genes that were previously shown to be consistently upregulated in mouse and human embryonic stem cells and embryonal carcinoma cells, as opposed to their differentiated counterparts ([@bib67]), and asked whether they are co-enriched in CTC clusters. First, a subset of 302 of these 335 genes were found to be expressed in our CTC samples from breast cancer patients (cutoff ≥ 3 transcripts per million). WGCNA with these genes identified four expression modules co-enriched in either single CTCs or CTC clusters. Particularly, we identified two co-expression modules, containing a total of 85 genes, enriched in patient-derived CTC clusters (n = 66 for blue; n = 19 for gray expression module) and in two co-expression modules, containing 217 genes, enriched in patient-derived single CTCs (n = 156 for turquoise; n = 61 for brown) ([Figures 4](#fig4){ref-type="fig"}A and 4B; [Table S4](#mmc4){ref-type="supplementary-material"}). GO network analysis of these modules confirmed that CTC clusters co-express stemness-related genes interconnected to cellular proliferation, while single CTCs co-express genes more related to metabolic processes ([Figures 4](#fig4){ref-type="fig"}C, 4D, and [S4](#figs4){ref-type="fig"}A; [Table S4](#mmc4){ref-type="supplementary-material"}). Importantly, TF target gene analysis ([@bib25]) of CTC cluster-enriched genes confirmed, among others, activity of SIN3A and OCT4, in line with our DNA methylation findings ([Figure 4](#fig4){ref-type="fig"}E; [Table S5](#mmc5){ref-type="supplementary-material"}). We then extended our analysis of the 302 stemness-related genes to xenograft-derived CTCs. In this case, WGCNA revealed four expression modules, of which two were found to be enriched in xenograft-derived CTC clusters and containing a total of 153 genes (n = 151 for green; n = 2 for yellow expression module) ([Figures 4](#fig4){ref-type="fig"}F and 4G; [Table S4](#mmc4){ref-type="supplementary-material"}), in high concordance (85% overlap) with patient-derived modules enriched in CTC clusters ([Figure 4](#fig4){ref-type="fig"}H; [Table S4](#mmc4){ref-type="supplementary-material"}). TF target gene analysis of those genes that are enriched in xenograft CTC clusters also revealed the activity of several stemness- and proliferation-related TFs, including SIN3A, OCT4, NANOG, and SOX2 ([Figure 4](#fig4){ref-type="fig"}I; [Table S5](#mmc5){ref-type="supplementary-material"}). Interestingly, the majority of these genes also displays hypomethylated promoter regions in CTC clusters compared to single CTCs ([Figures S4](#figs4){ref-type="fig"}B and S4C). Co-expression modules found in single CTCs from patients and xenografts also displayed significant overlap between each other (56.7%) ([Figure S4](#figs4){ref-type="fig"}D; [Table S4](#mmc4){ref-type="supplementary-material"}), and TF analysis revealed the activity of MYF6 and ASCL1, also displaying hypomethylated binding sites in single CTCs ([Figure S4](#figs4){ref-type="fig"}E; [Table S5](#mmc5){ref-type="supplementary-material"}).Figure 4Stemness-Related Gene Expression Analysis of Single CTCs and CTC Clusters from Breast Cancer Patients and Xenografts(A) WGCNA of 302 stemness-related transcripts in patient-derived single CTCs and CTC clusters, showing a hierarchical clustering tree of co-expression modules. Each module corresponds to a branch, which is labeled by a distinct color shown underneath.(B) WGCNA identifies four modules with highly correlated gene expression patterns in patient-derived single CTCs and CTC clusters. Correlations between each module and CTC clusters or single CTCs are indicated by the intensity of red or green color, respectively. p value for each module is shown in brackets.(C) GO network analysis of patient-derived transcripts identified in the CTC cluster-associated blue and gray modules using iRegulon. The node size indicates significance (p \< 0.05), and color intensity corresponds to the percentage of genes that are associated to each GO category. Indicative GO categories are shown.(D) GO term analysis of transcripts enriched in the blue and gray modules of patient-derived CTC clusters.(E) Patient-derived CTC cluster-associated blue and gray module gene regulatory network analysis showing putative transcription factor dependence on SIN3A, OCT4, and CBFB (green octagons).(F) WGCNA of 302 stemness-related transcripts in xenograft-derived single CTCs and CTC clusters, showing a hierarchical clustering tree of co-expression modules. Each module corresponds to a branch, labeled by a distinct color shown underneath.(G) WGCNA identifies four modules with highly correlated gene expression patterns in xenograft-derived single CTCs and CTC clusters. Correlations between each module and CTC clusters or single CTCs are indicated by the intensity of red or green color, respectively. p value for each module is shown in brackets.(H) Venn diagram showing the overlap of CTC cluster-enriched transcripts between patients (blue and gray modules) and xenograft (yellow and green modules).(I) Xenograft-derived CTC cluster-associated green and yellow module gene regulatory network analysis showing putative transcription factor dependence on SIN3A, OCT4, NANOG, BHLHE40, RORA, and FOXO1 (green octagons) and on SOX2 (orange circle).See also [Figure S4](#figs4){ref-type="fig"} and [Tables S4](#mmc4){ref-type="supplementary-material"} and [S5](#mmc5){ref-type="supplementary-material"}.Figure S4Stemness-Related Gene Expression Analysis of Single CTCs and CTC Clusters, Related to [Figure 4](#fig4){ref-type="fig"}(A) GO term analysis of transcripts identified in patient-derived single CTC-associated *turquoise* and *brown* modules.(B and C) Whisker plots showing the average methylation difference in CTC clusters relative to single CTCs, detected on the 5kb region upstream of the transcription start site of each target gene, in patient- (A) and xenograft-derived (B) CTCs. Hypomethylated genes in CTC clusters are represented with *blue* color, hypomethylated genes in single CTCs are represented with *red* color. Transcription factors relative to target genes are shown within *gray* boxes.(D) Venn diagram showing the overlap between genes enriched in single CTCs of patient-derived (*turquoise* and *brown*) and xenograft-derived (*orange* and *purple*) expression modules.(E) Regulatory network analysis of the 129 stemness-related genes that are commonly enriched in single CTCs in patient-derived (*turquoise* and *brown*) and xenograft-derived (*orange* and *purple*) expression modules, showing TF dependence on MYF6 and ASCL1, which also display hypomethylated binding sites based on the DMR analysis (*green* octagons).

Altogether, our gene expression data both at the transcriptome-wide level and also focused on stem cell-related genes strongly supports the model proposed with the DNA methylation analysis, suggesting that when compared to single CTCs, CTC clusters are endowed with a stemness- and proliferation-related network centered on the activity of key transcription factors including OCT4, SOX2, NANOG, and SIN3A. Activation of these programs may play a pivotal role in determining the metastasis-seeding ability of CTC clusters.

Identification of FDA-Approved Cluster-Targeting Agents {#sec2.5}
-------------------------------------------------------

Next, we sought to identify actionable vulnerabilities of CTC clusters, and to test whether the epigenetic and transcriptional features of clustered CTCs are reversible upon cluster dissociation into single cells. To this end, we evaluated 2,486 FDA-approved compounds ([Table S6](#mmc6){ref-type="supplementary-material"}) for their ability to dissociate clusters of human breast CTC-derived cells obtained from two patients with breast cancer (BR16 and BRx50), without affecting cellular viability. Cluster dissociation was assessed using an image-based high-content screening system and comparing CTCs treated with individual compounds to untreated CTCs as well as 40-μm-filtered single-cell suspension as negative and positive controls, respectively ([Figures 5](#fig5){ref-type="fig"}A, [S5](#figs5){ref-type="fig"}A, and S5B). For the majority of the 2,486 FDA-approved compounds, we observed no detectable reduction in cell viability (\>70% viability) or in the mean CTC cluster size (\> 450 μm^2^) of CTC-derived cells upon treatment (gray circles; [Figure 5](#fig5){ref-type="fig"}B). Yet we identified 39 compounds that significantly (p \< 0.0001) reduced mean CTC cluster size without compromising viability (orange circles; [Figure 5](#fig5){ref-type="fig"}B). These compounds include inhibitors of Na^+^/K^+^ ATPase (n = 6), histone deacetylase (HDAC) (n = 2), nucleotide biosynthesis (n = 5), kinases (n = 3), GPCRs (n = 2), cholesterol biosynthesis (n = 1), nuclear export (n = 1), tubulin (n = 10), DNA binding compounds (n = 8), and antibiotics (n = 1) ([Figures 5](#fig5){ref-type="fig"}B and [S5](#figs5){ref-type="fig"}C). Importantly, reducing compound concentration resulted in a gradual increase in mean CTC cluster size in both CTC-derived cell lines, followed by an increase in the number of nuclei detected, and a slight improvement in mitochondrial membrane potential and overall viability ([Figures 5](#fig5){ref-type="fig"}C and [S5](#figs5){ref-type="fig"}D). Under these conditions, six compounds consistently led to a significant decrease in mean CTC cluster size for both CTC-derived cell lines, even at the lowest concentrations tested. These compounds display similar functions and can be grouped into two families based on their mechanism of action, namely the Na^+^/K^+^-ATPase inhibitors digitoxin and ouabain octahydrate and the tubulin binding agents rigosertib, podofilox, colchicine, and vincristine sulfate ([Figures 5](#fig5){ref-type="fig"}C, [S5](#figs5){ref-type="fig"}C, and S5D).Figure 5Screen for FDA-Approved Compounds that Dissociate CTC Clusters(A) Representative images of unfiltered and filtered BR16 CTC-derived cells stained with Hoechst (blue) and TMRM (orange) (left). Representative images of single and clustered CTCs outline based on nuclei proximity as determined using the Colombus Image Analysis System (right). The plots show the mean CTC cluster size (area in micrometers squared) and percentage (%) of viability of unfiltered versus filtered BR16 cells (n = 4; ^∗∗∗^p \< 0.001 by Student's t test; ns, not significant) (bottom).(B) Effect of a 2-day treatment of BR16 cells with 2,486 FDA-approved compounds at 5 μM concentration, plotted as mean CTC cluster size (area in micrometers squared) versus percentage (%) of viability (n = 2). Thirty-nine FDA-approved compounds (orange circles) result in significant decrease in mean CTC cluster size (p \< 0.0001, F value = 7.71; DF = 38 using one-way ANOVA test; horizontal dashed red line, \< 450 μm^2^) and \> 70% detectable viability (vertical-dashed red line). BR16 cells that were untreated (red) or 40 μm filtered (green) are shown as controls.(C) The plot shows the mean CTC cluster size of BR16 cells treated with each of the 39 cluster-targeting compounds at four different concentrations. BR16 cells that were untreated (red) or 40 μm filtered (green) are shown as controls (top panel). The heatmap shows the number of nuclei, mean TMRM intensity, and percentage (%) of viability of BR16 cells treated with cluster-targeting compounds at the indicated concentrations (bottom). n = 2. Error bars represent SEM.See also [Figure S5](#figs5){ref-type="fig"} and [Table S6](#mmc6){ref-type="supplementary-material"}.Figure S5Screen for CTC Cluster-Targeting FDA-Approved Compounds, Related to [Figure 5](#fig5){ref-type="fig"}(A) Average TMRM intensity of unfiltered and filtered BR16 CTC-derived cells. n = 4; Error bars represent S.E.M; ^∗∗∗^p \< 0.001 by Student's t test.(B) Assessment of mean CTC cluster size (area in μm^2^), percent (%) of viability and TMRM intensity of unfiltered and filtered BRx50 CTC-derived cells. n = 4; Error bars represent SEM; ^∗∗∗^p \< 0.001 ^∗∗∗∗^p \< 0.0001 by Student's t test; ns, not significant.(C) Mean CTC cluster size and viability distribution of BR16 cells treated with cluster-targeting FDA-approved compounds at 5 μM concentration for 2 days (p \< 0.0001). BR16 cells that were not treated with FDA-approved compounds (*red*) or that were filtered (*green*) are shown for comparison. p = by one-way ANOVA test. n = 2; Error bars represent SEM.(D) Mean CTC cluster size (area in μm^2^) of BRx50 cells treated with 3 different concentrations (0.1, 0.5 and 1 μM) of CTC cluster-targeting FDA-approved compounds (n = 2; Error bars represent SEM). BRx50 cells that were not treated with FDA-approved compounds (*red*) or that were filtered (*green*) are shown for comparison (*top*). Heatmaps showing the number of nuclei, the average TMRM intensity and the percent (%) of viability of BRx50 cells treated with CTC cluster-targeting FDA-approved compounds at the indicated concentrations (*bottom*). Compounds that consistently dissociate CTC clusters in both BR16 and Brx-50 cells are highlighted in *red*.

Dissociation of CTC Clusters Leads to DNA Methylation Remodeling of Key Sites {#sec2.6}
-----------------------------------------------------------------------------

To assess whether CTC cluster dissociation into single cells would lead to DNA methylation remodeling, we treated four CTC-derived cell lines with the six CTC cluster-dissociating compounds individually. We found that a prolonged treatment (17 days) with 20 nM digitoxin or ouabain was unique in ensuring the highest degree of CTC cluster dissociation while allowing high viability and optimal proliferation rates comparable to untreated control cells ([Figure S6](#figs6){ref-type="fig"}A). Of note, the treatment period of 17 days was chosen to allow all cells (treated and untreated) to double at least 5 times and therefore enabling DNA methylation remodeling events to occur. Upon treatment, we processed BR16 and Brx50 cells for WGBS and RNA sequencing analysis to assess the molecular consequences of clusters dissociation. First, WGBS analysis revealed that the treatment did not affect global DNA methylation levels (data not shown). Yet a number of DMRs that were hypomethylated in CTC clusters gained methylation upon treatment ([Figures S6](#figs6){ref-type="fig"}B and S6C). Within the CTC cluster regions that gained methylation, i-*cis*Target analysis showed enrichment for binding sites of stemness-related TFs, such as OCT4, SOX2, and NANOG, and for SIN3A, among others ([Figure 6](#fig6){ref-type="fig"}A). Simultaneously, RNA expression analysis upon treatment with digitoxin or ouabain revealed a high concordance between the expression changes of stemness-related genes that were enriched in patient- and xenograft-derived CTC clusters (r = 0.69; p \< 2.22^−16^; [Figures 6](#fig6){ref-type="fig"}B), with the majority of these genes being downregulated as a consequence of the treatment (p \< 0.002; [Figure 6](#fig6){ref-type="fig"}B). TF analysis revealed that downregulated genes were targets of OCT4, NANOG, SOX2, and SIN3A transcription factors, among others ([Figure 6](#fig6){ref-type="fig"}C). In contrast, TF analysis of those genes that were upregulated as a consequence to CTC cluster dissociation highlighted target genes of other TFs, independent from a core pluripotency network ([@bib27]) ([Figure S6](#figs6){ref-type="fig"}D). Thus, CTC cluster dissociation into single cells with ouabain and digitoxin leads to DNA methylation remodeling and gain in methylation of critical binding sites for OCT4, SOX2, NANOG, and SIN3A, paralleled by downregulation of their corresponding target genes.Figure S6Treatment of CTC-Derived Cell Lines with Digitoxin and Ouabain, Related to [Figure 6](#fig6){ref-type="fig"}(A) Mean CTC cluster size (area in μm^2^) of BR16, BRx50, BRx07 and BRx68 CTC-derived cells treated with five different concentrations (1, 5, 10, 20 and 50 nM) of digitoxin or ouabain (n = 4; Error bars represent SEM). Cells that were not treated with FDA-approved compounds (*red*) or that were filtered (*green*) are shown for comparison (*top*). Heatmaps showing the number of nuclei, the average TMRM intensity and the percent (%) of viability of CTC-derived cells treated with digitoxin or ouabain at the indicated concentrations (*bottom*).(B and C) The histograms show the number of CTC cluster hypomethylated regions that become differentially methylated upon treatment of BR16 (B) and BRx50 (C) cells with 20 nM dixitoxin or ouabain (q ≤ 0.05). Regions that gain \>40% methylation are shown in *red*.(D) Regulatory network analysis of the genes (*orange*) that are commonly upregulated (log2 fold-change ≥ 0.5) in BR16 and BRx50 CTC-derived cell lines upon 17-day treatment with digitoxin or ouabain, showing TF dependence on PBX3, MAX, NFYA, ELF2, PARP1 and HOXA10. Upregulated genes that do not show dependence on the above TFs are shown below in *red*.(E) The plot shows the mean cell intensity of Ca^2+^-bound Fluo-3 after 30 min treatment of BR16 or BRx50 cells with 0.1, 0.5, 1 or 5 μM digitoxin or ouabain, respectively, relative to the untreated control (*red*) (n = 4; Error bars represent SEM).(F) Western blot for CLDN3, CLDN4 and GAPDH on BR16 cells with double knockout (2KO) of CLDN3 and CLDN4. The dashed line represents the point where irrelevant lanes were spliced out from the original scan.(G) Plot showing the reduction in mean cluster size (area in μm^2^) of the CLDN3/4 double KO BR16 cells, relative to control BR16 cells. ^∗^p \< 0.05; ^∗∗^p \< 0.01 by Student's t test. Error bars represent SEM.Figure 6Treatment with Na^+^/K^+^ ATPase Inhibitors and Tight Junction Dissociation Restores Methylation at Key Sites(A) i-*cis*Target analysis of CTC clusters hypomethylated regions, showing ≥ 40% methylation increase upon *in vitro* treatment with 20 nM of digitoxin or ouabain. Shown are enriched TFBSs in each CTC-derived cell line (NES ≥ 3).(B) The plot shows expression changes in stemness-related transcripts that are enriched in patient- and xenograft-derived CTC clusters (n = 168) upon treatment of BR16 and BRx50 cells with 20 nM digitoxin or ouabain. Pearson's correlation coefficient (r) and p values are shown.(C) Gene regulatory network analysis of downregulated genes upon treatment with 20 nM digitoxin or ouabain in both BR16 and BRx50 (log FC ≤ −0.5), showing putative dependence on OCT4, SOX2, NANOG, and SIN3A, among other TFs (green octagons). Downregulated genes that are not regulated by the TFs mentioned above are shown in red.(D) The plot shows the mean CTC cluster size of BR16 and BRx50 CTC-derived cells treated with FCCP or CCCP at different concentrations (n = 4). BR16 and BRx50 cells that were untreated (red) or 40 μm filtered (green) are shown as controls (top). The heatmap shows the percentage (%) of viability of treated cells at the indicated concentrations (middle). The plot shows the mean cell intensity of Ca^2+^-bound Fluo-3 after treatment with FCCP or CCCP, relative to the untreated control (red) (n = 4) (bottom). Error bars represent SEM.(E) Heatmap showing methylation variable regions among CTC cluster hypomethylated regions that gain ≥ 40% methylation upon CLDN3/4 double knockout in BR16 cells.(F) i-*cis*Target analysis of CTC cluster hypomethylated regions, showing ≥ 40% methylation increase upon CLDN3/4 double knockout in BR16 cells. Shown are enriched TFBSs (NES ≥ 3).See also [Figure S6](#figs6){ref-type="fig"} and [Tables S2](#mmc2){ref-type="supplementary-material"}, [S3](#mmc3){ref-type="supplementary-material"}, and [S7](#mmc7){ref-type="supplementary-material"}.

Inhibition of the Na^+^/K^+^ ATPase using ouabain and intracellular increase of calcium levels have been previously shown to negatively affect the formation of tight junctions and desmosomes in epithelial cells ([@bib50], [@bib54], [@bib61]). Here, to gain more insights into the mechanism of action of Na^+^/K^+^ ATPase inhibitors in the context of CTC clusters disruption, we tested whether the suppression of ATP production and consequently, an increase in intracellular calcium levels---using the proton uncouplers FCCP and CCCP ([@bib12])---would also dissociate CTC clusters. Our results show that treatment of CTC-derived cells with increasing concentrations of FCCP and CCCP for 18 hr leads to a gradual decrease of mean cluster size ([Figure 6](#fig6){ref-type="fig"}D), paralleled by a significant increase in intracellular Ca^2+^ levels ([Figure 6](#fig6){ref-type="fig"}D). Interestingly, a similar gradual increase in intracellular Ca^2+^ is also observed when treating with increasing concentrations of digitoxin and ouabain ([Figure S6](#figs6){ref-type="fig"}E). These results support a model whereby intracellular Ca^2+^ increase leads to the inability of cancer cells to properly form cell-cell junctions.

To further confirm this model, we assessed whether cell-cell junction disruption in CTC-derived cells would lead to clusters dissociation as well as DNA methylation remodeling at CTC cluster-associated DMRs. To this end, we employed the CRISPR technology to simultaneously knockout both claudin 3 (CLDN3) and claudin 4 (CLDN4) in BR16 CTC-derived cells, two of the highest-expressed tight junction proteins in CTC clusters ([Table S4](#mmc4){ref-type="supplementary-material"}). Using two independent sgRNAs for each gene, we generated three BR16 lines with double CLDN3/4 knockout, which also displayed a significant reduction of mean CTC cluster size ([Figures S6](#figs6){ref-type="fig"}F and S6G). Whole-genome bisulfite sequencing of the CLDN3/4 double-knockout cells showed that, upon dissociation into single cells and similarly to the events that occurred upon Na^+^/K^+^ ATPase inhibition, a number of CTC cluster-associated hypomethylated regions gained methylation ([Figure 6](#fig6){ref-type="fig"}E). Interestingly, i-*cis*Target analysis of the regions that gained higher levels of methylation revealed an enrichment of binding sites for OCT4, SOX2, NANOG, and SIN3A ([Figure 6](#fig6){ref-type="fig"}F), further indicating that CTC clustering directly affects DNA methylation dynamics at bindings sites for stemness- and proliferation-associated TFs.

Together, our results indicate that Na^+^/K^+^ ATPase inhibition leads to CTC clusters dissociation through the increase of the intracellular Ca^2+^ concentration and the consequent inhibition of cell-cell junction formation, resulting in DNA methylation remodeling at critical stemness- and proliferation-related binding sites.

Treatment with Na^+^/K^+^-ATPase Inhibitors Suppresses Spontaneous Metastasis Formation {#sec2.7}
---------------------------------------------------------------------------------------

To test whether ouabain and digitoxin would also enable CTC clusters disruption *in vivo*, we took a dual approach. First, we tested whether a 17-day *in vitro* treatment with ouabain and digitoxin would translate into a reduced ability of the treated cells to efficiently seed metastasis in untreated mice ([Figure 7](#fig7){ref-type="fig"}A). To this end, upon treatment, BR16 cells stably expressing GFP-luciferase were injected into the tail vein of NSG mice and noninvasively monitored through luminescence imaging for their ability to seed and propagate metastatic lesions. We found that while the treatment with digitoxin or ouabain did not affect the ability of BR16 cells to lodge in the lung tissue immediately after injection (see "day 0"; [Figure 7](#fig7){ref-type="fig"}B), it led to a reduced ability to survive during the first day upon arrival, as confirmed by a significant increase in the expression of cleaved caspase 3 compared to control cells (see "day 1"; [Figures 7](#fig7){ref-type="fig"}B, [S7](#figs7){ref-type="fig"}A, and S7B). Overall, this difference in the ability to survive during the very early steps of metastasis seeding resulted in a delayed metastatic outgrowth despite the absence of further treatment *in vivo*, as measured during the course of 72 days upon injection ([Figure 7](#fig7){ref-type="fig"}C).Figure 7Treatment with Na^+^/K^+^ ATPase Inhibitors Suppresses Spontaneous Metastasis Formation(A) Schematic representation of the experiment.(B) The plots show the total bioluminescence flux at day 0 (left) and day 1 (right) upon tail vein injection of BR16 cells pre-treated with 20 nM digitoxin or ouabain. n = 5 for controls and ouabain, n = 4 for digitoxin; ^∗^p \< 0.05 by Student's t test. ns, not significant. Error bars represent SEM.(C) Metastasis growth curve over 72 days upon tail vein injection of BR16 cells pre-treated with 20 nM digitoxin or ouabain. n = 5; ^∗^p \< 0.05; ^∗∗^p \< 0.01 by Student's t test. Error bars represent SEM.(D) Schematic representation of the experiment.(E) The plots show the percentage (%) of spontaneously generated single CTCs and CTC clusters detected in the blood of BR16 xenografts treated with ouabain. n = 5; ^∗∗∗^p \< 0.001 by Student's t test. error bars represent SEM.(F) The plot shows the metastatic index of BR16 xenografts treated with ouabain. n = 11 for controls, n = 5 for ouabain; ^∗∗^p \< 0.01 by Student's t test. Error bars represent SEM.(G) Representative images of the bioluminescence signal measured in brain and in liver of control and ouabain-treated NSG mice.See also [Figure S7](#figs7){ref-type="fig"}.Figure S7Treatment with Digitoxin and Ouabain Reduces Metastasis Formation, Related to [Figure 7](#fig7){ref-type="fig"}(A) The plots show the percent of Ki67-positive cancer cells detected in the lungs of NSG mice at Day 0 or Day 1 upon injection with BR16 CTC-derived cells, treated *in vitro* with digitoxin or ouabain. Cancer cells are identified through Pan Cytokeratin staining; n = 4 mice for each condition. Error bars represent SEM; ns, not significant.(B) The plots show the percent of Caspase 3-positive cancer cells detected in the lungs of NSG mice at Day 0 or Day 1 upon injection with BR16 CTC-derived cells, treated *in vitro* with digitoxin or ouabain. Cancer cells are identified through Pan Cytokeratin staining; n = 4 mice for each condition. ^∗^p \< 0.05 by Student's t test; Error bars represent SEM; ns, not significant.(C) The plot shows the total bioluminescence flux emitted from the primary tumor of BR16 xenografts treated with vehicle (control) or ouabain. Error bars represent SEM; ns, not significant.(D) The plot shows the total number of CTCs, including both single CTCs and CTC clusters, detected per mL of blood in BR16 xenografts treated with vehicle(control) or ouabain. n = 6 for controls and n = 5 ouabain; Error bars represent SEM; ns, not significant.(E) The plots show the percent (%) of spontaneously generated single CTCs and CTC clusters detected in the blood of LM2 xenografts treated with vehicle (control) or ouabain. n = 11 for controls, n = 8 for ouabain; ^∗∗^p \< 0.01.(F) The plot shows the total bioluminescence flux emitted from the primary tumor of LM2 xenografts treated with vehicle (control) or ouabain. Error bars represent SEM; ns, not significant.(G) The plot shows the total number of CTCs, including both single CTCs and CTC clusters, detected per mL of blood in LM2 xenografts treated with vehicle (control) or ouabain. n = 11 for controls and n = 8 ouabain; Error bars represent SEM; ns, not significant.(H) The plot shows the metastatic index of LM2 xenografts treated with vehicle (control) or ouabain. n = 18 for controls, n = 8 for oubain. ^∗^p \< 0.01 by Student's t test; Error bars represent SEM.

Second, mimicking more closely the clinical setting, to assess the effect of our CTC cluster-dissociation strategy for the spontaneous formation of CTC clusters and metastasis from a primary tumor, we injected BR16 cells in the mammary fat pad of NSG mice. Fourteen weeks after primary tumor formation we administered ouabain daily for three weeks and assessed CTC composition and the occurrence of spontaneous metastatic lesions ([Figure 7](#fig7){ref-type="fig"}D). Importantly, we observed that ouabain treatment reduced the frequency of spontaneously generated CTC clusters while increasing the frequency of single CTCs ([Figure 7](#fig7){ref-type="fig"}E), without altering the size of the primary tumor nor overall CTC numbers ([Figures S7](#figs7){ref-type="fig"}C and S7D). Along with a reduction in the frequency of CTC clusters, ouabain treatment also resulted in a remarkable suppression (80.7-fold) of the total metastatic burden ([Figures 7](#fig7){ref-type="fig"}F and 7G). In a similar fashion, when administering ouabain treatment to NSG mice carrying spontaneously metastasizing LM2 tumors, we also observed an increase in the proportion of single CTCs and a decrease in CTC clusters ([Figure S7](#figs7){ref-type="fig"}E), without any change in the primary tumor size nor overall CTC numbers ([Figures S7](#figs7){ref-type="fig"}F and S7G), leading to a reduced metastatic burden compared to control ([Figure S7](#figs7){ref-type="fig"}H).

Together, these results demonstrate that Na^+^/K^+^ ATPase inhibition *in vivo* suppresses the ability of a cancerous lesion to spontaneously shed CTC clusters, leading to a remarkable reduction in metastasis seeding ability.

Discussion {#sec3}
==========

Our study provides a comprehensive genome-wide analysis of the DNA-methylation events that characterize CTCs in patients and xenografts. Surprisingly, we find that phenotypic differences---such as the ability of CTCs to navigate through the bloodstream as single cells or multicellular clusters---shape the DNA methylome. Clustering of CTCs results in hypomethylation of binding sites that are typically occupied by master stemness and proliferation regulators, including OCT4, NANOG, SOX2, and SIN3A, and hypermethylation of Polycomb target genes. More globally, we also find that the DNA methylation profile of CTC clusters is detected at the level of the primary tumor in a subset of breast cancers that are characterized by a poor prognosis. CTC clusters dissociation into single cells with Na^+^/K^+^ ATPase inhibitors or through cell-cell junction knockdown enables DNA methylation remodeling at critical sites, highlighting a direct connection between clustering and methylation status. As a result, Na^+^/K^+^ ATPase inhibitors treatment emerges as a new strategy to significantly reduce the spread of cancer, providing a rationale for using these compounds in clinical studies.

Our results suggest that CTC clusters may share several properties that commonly feature stem cell biology. For instance, OCT4, NANOG, SOX2, and SIN3A are predominantly active in embryonic stem cells (ESCs), simultaneously regulating self-renewal and proliferation ([@bib27], [@bib47], [@bib51], [@bib62]). In addition, ESCs rely on Polycomb-mediated repression of differentiation genes and chromatin remodeling to maintain their active pluripotency network ([@bib36]). Cell-cell junction activity has been shown in several instances to safeguard pluripotency and to be required for a complete reprogramming of somatic cells into stem cells, and disruption of cell-cell junctions (e.g., through targeting of E-cadherin) in human ESC results into OCT4, NANOG, and SOX2 downregulation along the loss of stemness features ([@bib38], [@bib39], [@bib40], [@bib53]). Thus, by analogy with stem cell biology, elevated expression of cell-cell junction components in cancer cells may not only enable their intravasation in the bloodstream as multicellular clusters but also their ability to retain stem-like features that facilitate metastasis initiation.

Identifying FDA-approved compounds that dissociate CTC clusters provides a new tool to reduce the spread of cancer. Our results demonstrate that CTC clusters disruption into single cells via the inhibition of the Na^+^/K^+^ ATPase has a dual, yet interconnected effect. The fact that clusters dissociation leads to DNA methylation remodeling at critical sites provides a direct link between a phenotypic state of the cells (i.e., clustered versus single) and DNA methylation dynamics. Second, CTC clusters disruption increases the proportion of single CTCs in the bloodstream but suppresses overall metastasis formation, indicating that targeting CTC clusters could be a valuable therapeutic strategy. In this regard, while we cannot exclude that an anti-cluster treatment might be beneficial also at the late disease stages in those patients whose metastasis may seed other metastases ([@bib48], [@bib57], [@bib70]), *in vivo* data suggest that treatment with ouabain and digitoxin should be administered early, ideally at the time of localized disease and before dissemination to distant sites, with the objective to prevent CTC cluster formation.

Ouabain and digitoxin are cardiac glycosides with a similar chemical structure, used in low doses for the treatment of hypotension and cardiac arrhythmias, acting through nonselective inhibition of the Na^+^/K^+^ ATPase ([@bib7], [@bib59]). Upon inhibition of the Na^+^/K^+^ ATPase, cellular uptake of Na^+^ occurs, leading to a simultaneous increase in intracellular Ca^2+^ and impaired translocation of desmosomal and tight junction proteins to the cellular membrane ([@bib7], [@bib8]). Based on this, our experimental evidences support a model whereby pharmacological inhibition of the Na^+^/K^+^ ATPase either directly or indirectly (i.e., by depletion of available ATP) leads to a concomitant increase in intracellular Ca^2+^ levels and CTC cluster disruption through suppression of functional cell-cell junction assembly ([@bib10], [@bib29]).

Together, our study provides key insights into the biology of CTCs and highlights a fundamental connection between phenotypic features of CTCs (such as their ability to circulate as multicellular clusters) and DNA methylation dynamics at critical stemness- and proliferation-related sites. Further, we identify the Na^+^/K^+^ ATPase inhibitors ouabain and digitoxin as FDA-approved agents capable to dissociate CTC clusters into single cells and to suppress spontaneous metastasis formation in xenograft models, providing a rationale for applying these compounds for the treatment of patients with breast cancer.

STAR★Methods {#sec4}
============

Key Resources Table {#sec4.1}
-------------------

REAGENT or RESOURCESOURCEIDENTIFIER**Antibodies**EpCAM-AF488Cell SignalingCat\# CST5198; RRID: [AB_10692105](nif-antibody:AB_10692105){#interref85}HER2-AF488BioLegendCat\# 324410; RRID: [AB_2099256](nif-antibody:AB_2099256){#interref90}EGFR-FITCGeneTexCat\# GTX11400; RRID: [AB_368217](nif-antibody:AB_368217){#interref95}CD45-BV605 (anti human)BiolegendCat\# 304042; RRID: [AB_2562106](nif-antibody:AB_2562106){#interref100}CD45-BV605 (anti mouse)BiolegendCat\# 103140; RRID: [AB_2562342](nif-antibody:AB_2562342){#interref105}**Chemicals, Peptides, and Recombinant Proteins**2,486 FDA-approved moleculesNEXUS platform, ETH ZurichRefer to [Table S6](#mmc6){ref-type="supplementary-material"} for individual chemical abstracts service (CAS) registry \#**Deposited Data**RNA sequencingThis paperGEO: [GSE111065](ncbi-geo:GSE111065){#intref0010}Bisulfite sequencingThis paperENA: PRJEB25101**Experimental Models: Cell Lines**BR16 CTC-derivedAceto lab, University of Baseln/aBrx50, Brx07 and Brx68 CTC-derivedHaber lab, MGH Cancer Centern/aMDA-MB-231 LM2Massague lab, MSKCCn/a

Contact for Reagent and Resource Sharing {#sec4.2}
----------------------------------------

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Nicola Aceto (<Nicola.Aceto@unibas.ch>).

Experimental Model and Subject Details {#sec4.3}
--------------------------------------

### Human blood samples collection {#sec4.3.1}

Blood specimens for CTC analysis were obtained from the University Hospital Basel after informed patient consent, according to protocols EKNZ BASEC 2016-00067 and EK 321/10, which received ethical approval from the Swiss authorities (EKNZ, Ethics Committee northwest/central Switzerland).

### Mouse blood samples collection {#sec4.3.2}

For mouse studies, blood was retrieved via cardiac puncture and up to 1 mL of blood was collected. All mouse experiments were carried out in compliance with institutional and cantonal guidelines (approved mouse protocol \#2781, cantonal veterinary office of Basel-City).

### Cell lines {#sec4.3.3}

BR16 CTC-derived cells were generated from a patient with metastatic breast cancer at the University Hospital Basel as previously described ([@bib69]). Brx50, Brx07 and Brx68 CTC-derived cells were obtained from the Haber and Maheswaran lab (MGH Cancer Center, Harvard Medical School, Boston, MA). MDA-MB-231 (LM2) cells were donated from Joan Massague's lab (MSKCC, New York, NY, USA).

Method Details {#sec4.4}
--------------

### Cell Culture {#sec4.4.1}

CTC-derived cells were maintained under hypoxic conditions (5% oxygen) on ultra-low attachment (ULA) 6-well plates (Corning, Cat\# 3471-COR). CTC growth medium containing 20 ng/ml recombinant human Epidermal Growth Factor (GIBCO, Cat\# PHG0313), 20 ng/ml recombinant human Fibroblast Growth Factor (GIBCO, Cat\#100-18B), 1x B27 supplement (Invitrogen, Cat\#17504-044) and 1x Antibiotic-Antimycotic (Invitrogen, Cat\# 15240062) in RPMI 1640 Medium (Invitrogen, Cat\# 52400-025) was added every third day. For passaging, cells were spun down at 800 g for 5 min using a Heraeus Multifuge X3R centrifuge (Invitrogen, Cat\#75004515). The supernatant was subsequently aspirated and cells were resuspended in 2 ml/well CTC medium and plated in 6-well ULA plates. LM2 cells were passaged in DMEM/F-12 medium (Invitrogen, Cat\#11330057) supplemented with 10% FBS (Invitrogen, Cat\# 10500064) and Antibiotic-Antimycotic (Invitrogen, Cat\# 15240062). For passaging, LM2 cells were washed once with D-PBS (Invitrogen, Cat\#14190169) and dissociated using 0.25% Trypsin (Invitrogen, Cat\#25200056).

### CTC Capture and Identification {#sec4.4.2}

An average of 7.5 mL of blood per patient or up to 1 mL per mouse was drawn in EDTA vacutainers. Within 1 hr from blood draw, the blood was processed through a Parsortix GEN3D6.5 Cell Separation Cassette (Angle Europe) according to the manufacturer's instructions. Captured CTCs were further stained on the Parsortix cassette with EpCAM-AF488 (CellSignaling, Cat\# CST5198), HER2-AF488 (\#324410, BioLegend), EGFR-FITC (GeneTex, Cat\# GTX11400) and CD45-BV605 (Biolegend, Cat\# 304042 (anti-human); Cat\# 103140 (anti-mouse)) antibodies or imaged directly for GFP.

### Tumorigenesis Assays {#sec4.4.3}

For tail vein experiments, NOD SCID Gamma (NSG) mice (Jackson Labs) were injected with 7x10^5^ BR16-mCherry cells resuspended in 100 μL D-PBS and monitored with IVIS Lumina II (Perkin Elmer). For xenografts, 1x10^6^ LM2-GFP, 1x10^6^ BRx50-GFP or 1x10^6^ BR16-GFP cells were resuspended in 100 μL of 50% Cultrex PathClear Reduced Growth Factor Basement Membrane Extract (R&D Biosystems, Cat\# 3533-010-02) in D-PBS and injected orthotopically in NSG mice. Blood draw was performed 4-5 weeks after tumor onset for LM2 cells, 5 months after tumor onset for BR16 and 6-7 months after tumor onset for BRx50 cells. For *in vivo* treatment with ouabain, NSG mice were injected IP daily with 100 μL containing 0.67 mg/Kg ouabain diluted in PBS.

### Immunofluorescence Staining {#sec4.4.4}

For immunofluorescence staining of CTCs, after Parsortix enrichment, CTCs were released in PBS and cytospinned for 3 min at 500 rpm on coated microscope slides (Thermo Scientific, Cat\# 5991056), left to air dry for 5 min and then fixed with 4% PFA (Sigma; Cat\# 58127) for 12 min. For lung tissue, immunoflourescence was performed on 4 μm paraffin sections. Paraffin sections were deparaffinized prior to immersion in Tris-EDTA buffer (10 mM Tris base/1 mM EDTA/0.05% Tween20, pH 9.0) at 95°C for 40 minutes for antigen retrieval. Slides were then washed twice in PBS and premeabilized with 5% Triton X-100 (Sigma) in PBS for 20 min, washed once with 0.2% Tween in PBS and twice with PBS. Blocking was performed for 1 hour at room temperature in 10% horse serum (Sigma, Cat\# 2-05F00-I) followed by 2 times wash in PBS. Primary antibodies for Pan CytoKeratins (Abcam, Cat\# ab190707), Ki67 (Abcam; Cat\# ab15580) and Cleaved Caspase-3 (Cell Signaling, Cat\# 9664S) were added at 1:50, 1:60 and 1:100 dilutions, respectively, in 0.5% BSA in PBS and incubated overnight at 4°C. The following day the slides were washed once with 5% horse serum, 0.2% tween in PBS followed by one wash with PBS. Secondary antibodies anti mouse-AF488 (Invitrogen; Cat\# A21202) and anti-rabbit-AF647 (Invitrogen; Cat\# A31573) were added at 1:200 dilution in 0.5% BSA in PBS for 30 min followed by three washes with PBS. Finally, slides were incubated for 5 min with 10 μM DAPI in PBS (Sigma, Cat\# D9542) and mounted using ProLong Gold Antifade (Invitrogen; Cat\# P36934).

### FCCP, CCCP treatment and Intracellular Ca^2+^ measurement {#sec4.4.5}

FCCP (Abcam, Cat\# ab120081) and CCCP (Abcam, Cat\# ab141229) were reconstituted in DMSO at 100 mM concentration and further diluted in CTC medium to the indicated concentrations. For Ca^2+^ measurements, Fluo-3 AM (cell permeant, Invitrogen; Cat\# F1241) was added to the 20'000 cells/well into a 96 well Black/clear Tissue culture treated plate (BD Falcon, Cat\#353219) at 2 μM concentration together with 4 μM Hoechst 34580 (Invitrogen, Cat\# H21486) and 4 μM TOTO-3 (Invitrogen, Cat\# T3604) and cells were incubated for 30 min at 37°C. Plates were scanned using Operetta High Content Imaging System (Perkin Elmer) and Ca^2+^ bound Fluo-3 intensity per cell was calculated using Columbus Image Data Storage and Analysis System (Perkin Elmer).

### CRISPR-CAS9 CLDN3/4 Double knock out in BR16 {#sec4.4.6}

We used lentiviral delivery of pLenti-Cas9-EGFP vector (Addgene) to generate a BR16 CTC-derived cell line that stably expresses the Cas9 protein together with GFP. In BR16-Cas9-GFP line we then introduced sgRNA sequences that target either CLDN3 or CLDN4. In detail, sgRNA sequences were designed using the GPP Web Portal (<https://portals.broadinstitute.org/gpp/public/analysis-tools/sgrna-design>). Two sgRNAs targeting CLDN3 ((sense) 5′-CACGTCGCAGAACATCTGGG-3′ and (sense) 5′-ACGTCGCAGAACATCTGGGA-3′) were cloned in vector pLentiGuide-Puro (Addgene) and 2 sgRNAs targeting CLDN4 ((sense) 5′-CAAGGCCAAGACCATGATCG-3′ and (sense) 5′-ATGGGTGCCTCGCTCTACGT-3′) were cloned in vector pLentiGuide-Blast. Vector pLentiGuide-Blast was generated by replacing puromycin resistance gene on plasmid pLentiGuide-Puro with the blasticidin resistance gene using the MluI and BsiWI restriction enzyme sites. Double positive-clones were selected based on puromycin (1 μg/ml) and blasticidin (10 μg/ml) antibiotic selection for 2 weeks and CLDN3/CLDN4 knockout was verified by western blot.

### Western Blot {#sec4.4.7}

Cells were pelleted at 800 rpm for 5 min and lysed in 1% SDS, 50 mM Tris pH7.8, 10 mM EDTA pH 8.0 and 0.1 M DTT in PBS with the addition of protease inhibitors according to manufacturer's instructions (Complete EDTA-free Protease Inhibitor Cocktail Tablets - Roche; Cat\#11873580001). Cell lysates were passed through a 27-G syringe to achieve homogeneous disruption and centrifuged at 10'000 rpm for 5 min at 4°C. The protein containing supernatant was transferred to a fresh tube. Protein concentration was determined using the DC Protein Assay Detection System (Biorad, Cat\# 500-0114). The proteins were separated on a 12% SDS gel using the Bolt system (Invitrogen) and transferred to a nitrocellulose membrane using the iBlot2 system (Invitrogen) according to manufacturer's instructions. Membranes were blocked in 5% milk in PBS for 1 hour at room temperature and primary antibodies for the detection of CLDN3 (1:1000 dilution; Novus; Cat\#NB100-91711), CLDN4 (1:200 dilution; Invitrogen; Cat\# 32-9400) and GAPDH (1:1000 dilution; Cell Signaling; Cat\# 2118S) were added for incubation overnight at 4°C. The next day membranes were washed 3 times for 5 min in 0.2% tween in PBS. Secondary antibodies were added at 1:5000 dilution in 5% Milk in PBS and detection was performed using the Odyssey imaging system (Li-Cor) all according to manufacturer's instructions.

### Single-Cell Micromanipulation {#sec4.4.8}

CTCs were harvested from the Parsortix cassette in 1 mL D-PBS solution (Invitrogen, Cat\#14190169) in a 6-well ultra-low attachment plate (Corning, Cat\# 3471-COR) and visualized using a CKX41 Olympus inverted fluorescent microscope (part of the AVISO CellCelector Micromanipulator --ALS). Single CTCs and CTC clusters were identified based on intact cellular morphology, AF488/FITC-positive staining and lack of BV605 staining. Target cells were individually micromanipulated with a 30 μm glass capillary on the AVISO CellCelector micromanipulator (ALS) and deposited into individual PCR tubes (Axygen, Cat\#321-032-501) containing 10 μL of 2x Digestion Buffer (EZ DNA Methylation Direct Kit -Zymo, Cat\# D5020) for WGBS or 2 μL of RLT lysis buffer (QIAGEN, Cat\#79216) supplemented with 1U/μl SUPERase In RNase inhibitor (Invitrogen, Cat\# AM2694) for RNA sequencing, and immediately flash frozen in liquid nitrogen.

### Single Cell Whole-genome Bisulfite Sequencing {#sec4.4.9}

Proteinase K digestion and bisulfite treatment was performed according to the manufacturer's instructions for EZ DNA Methylation Direct Kit (Zymo, Cat\# D5020). Bisulfite-treated DNA was eluted using 9 μL of Elution Buffer and used for library generation with TruSeq DNA methylation kit (Illumina, Cat\# EGMK91396) according to manufacturer's instructions. For amplification, 18 cycles were performed using Failsafe Enzyme (Illumina, Cat\# FSE51100) and indexes were introduced with Index Primers' Kit (Illumina, Cat\# EGIDX81312). Library purification was performed using Agencourt AMPure XP beads at a ratio of 1:1 according to manufacturer's instructions. To avoid DNA loss during pipetting steps, Corning DeckWork low binding barrier pipet tips were used (Sigma, Cat\# CLS4135-4X960EA). Library concentration was estimated using Qubit DS DNA HS Assay Kit according to manufacturer's instructions (Invitrogen, Cat\#Q32854).

### RNA-Seq Library generation {#sec4.4.10}

RNA was captured on beads conjugated with oligo-dT primer according to a published protocol ([@bib45]). cDNA was generated according to the Smart-Seq 2 protocol ([@bib52]). Sequencing libraries were generated and indexed from 0.25 ng of cDNA per sample using the Nextera XT DNA Library Preparation Kit (Illumina, Cat\# FC-131-2001) according to manufacturer's instructions.

### FDA-Approved Compound Screen {#sec4.4.11}

A library containing 2,486 FDA-approved compounds was obtained from the Nexus Platform -- ETH Zurich. Each compound was resuspended using CTC medium at a 15 μM concentration and 20 μL were aliquoted in duplicate in a total of 64 Flat Bottom Clear Ultra Low attachment 96-well plates (Corning, Cat\#3474). To obtain a single cell suspension from CTC-derived cells, a 40 μm cell strainer was used (Corning, Cat\# 431750). 40 μL containing 15,000 CTC-derived cells were seeded per well in 96-well ultra-low attachment plates that contained 20 μL of pre-aliquoted FDA-approved compounds at 15 μM concentration, so that the final compound concentration was 5 μM. Plates were incubated in hypoxic conditions (5% oxygen) for 2 days and then 20 μL were transferred into a 96 well Black/clear Tissue culture treated plate (BD Falcon, Cat\#353219) containing 40 μL of D-PBS (Invitrogen, Cat\#14190169) and stained for 1hr at 37°C with a final concentration of 4 μM Hoechst 34580 (Invitrogen, Cat\# H21486), 2 μM TMRM (Invitrogen, Cat\#T668) and 4 μM TOTO-3 (Invitrogen, Cat\# T3604). For each plate, two negative controls (non-treated cells) and two positive controls (40 μm-filtered cells) were included. Z-factors were calculated per individual plate using the following formula: Z′ = 1-3(σ~s~ + σ~c~)/\|μ~s~-μ~c~\| (σ: standard deviation, μ: mean, s: positive control and c: negative control) ([@bib46]) and ranged between 0.62-0.937. Plates were scanned using Operetta High Content Imaging System (Perkin Elmer) and CTC cluster analysis was performed using Columbus Image Data Storage and Analysis System (Perkin Elmer).

Quantification and Statistical Analysis {#sec4.5}
---------------------------------------

### Whole Genome Bisulfite Data Processing {#sec4.5.1}

After sequencing, initial quality assessment of data was performed using FastQC (<https://www.bioinformatics.babraham.ac.uk/projects/fastqc>). Adaptor sequences, low quality ends, first 9 base pairs from 5′ end and 3 base pairs from 3′ were removed with Trim Galore (v0.4.2, <http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/>; parameters:--paired--retain_unpaired -r1 36 -r2 36--clip_r1 9--three_prime_clip_r1 3--clip_r2 9--three_prime_clip_r2 3--quality 20--phred33). Trimmed reads were aligned to GRCh38 human genome using Bismark Bisulfite Mapper (v0.17.0; paramters:--non_directional--bowtie2) ([@bib32]) with the alignment tool Bowtie2 (v2.2.9) ([@bib33]). Reads with a conversion rate lower than 90% or less than three non-converted cytosines in non-CpG context were removed. Finally, methylation calls were extracted after deduplication using Bismark. Resulting data was processed using MethylKit ([@bib3]) (v1.4.0) to aggregate methylation values into genomic tiling regions of 5 kilobases using a step size of 1kilobase normalizing coverage between samples. Annotation used for metaplots of genes and CpG islands was performed using genomation toolkit in R/bioconductor package ([@bib4]).

### Differential Methylation Analysis and Enrichment Analysis {#sec4.5.2}

Differential methylation at genomic tiles was assessed using Methylkit Chisq method to compare the fraction of methylated CpGs in test and control conditions, and correcting for overdispersion. Tiles with at least 2 samples covered per group and with a Q-value ≤ 0.05 using the SLIM ([@bib64]) method were considered as differentially methylated. The region set enrichment analysis was performed following the example of [@bib18]. Normalized enrichment scores (NESs) of TBFSs in DMRs were calculated using i-*cis*Target (<https://gbiomed.kuleuven.be/apps/lcb/i-cisTarget/>) ([@bib22]) using full motif analysis and NES score threshold ≥ 3. Gene ontology network analysis was generated using ClueGO ([@bib63]). LOLA (v1.8; <http://databio.org/regiondb>) ([@bib60]) was used to annotate the methylation data with a catalog of regulatory region sets from Codex database and ENCODE.

### Enrichment Analysis Based on Genomic Features {#sec4.5.3}

Individual CpGs were mapped to genes and their promoters using RefSeq genes annotation from the UCSC genome browser (<https://genome.ucsc.edu>; date: 01/06/2018). Promoters were defined as ± 2 kb region around the transcription start site (TSS). Mapping to super-enhancer regions was based on dbSUPER (<http://asntech.org/dbsuper/>), an integrated database of super-enhancers that provides a list of genes associated to each region. Genomic coordinates from dbSUPER were converted from GRCh37 to GRCh38 via LiftOver (<https://genome.ucsc.edu/cgi-bin/hgLiftOver>). Each genomic feature was interrogated for differential methylation in the same fashion as we did for genomic tiles. Genes corresponding to genomic features with a nominal *P value* ≤ 0.01 and an absolute methylation difference ≥ 20% were considered for enrichment analysis using the gene set collection C2 from the Molecular Signatures Database v6.1 (MSigDB). Gene set enrichment analysis was performed using hypergeometric test implemented in the clusterProfiler R/bioconductor package. Gene-sets with an adjusted *P value* ≤ 0.05 were considered significant.

### Survival Analysis using TCGA DNA Methylation Data {#sec4.5.4}

A total of 789 primary breast cancer samples from TCGA were available from the Illumina HumanMethylation 450 (HM450K) array. Level 3 data was downloaded from The National Cancer Institute (NCI) Genomic Data Commons (GDC) using the TCGAbiolinks bioconductor package ([@bib11]). Downloaded data are β-values ranging from 0 to 1 and reflecting the fraction of methylated alleles at each CpG. Probes overlapping regions of the CTC cluster-associated DNA methylation signature were selected for further analysis. A curated set of clinical parameters was obtained from the TCGA Pan-Cancer Clinical Data Resource (TCGA-CDR) ([@bib43]). For survival analysis, participants were divided into quartiles based on the average methylation across the CTC cluster-associated DNA methylation signature loci and progression free survival (PFS) was used as an endpoint. The Kaplan-Meier survival curves were drawn comparing the lowest (25%, Q1) versus the highest (75%, Q4) quartile and log-rank test *P value* was calculated using the survival R package.

### RNA-Seq Data Processing {#sec4.5.5}

Initial quality assessment for RNA-seq data was performed using FastQC (<https://www.bioinformatics.babraham.ac.uk/projects/fastqc>), FastQ Screen (<https://www.bioinformatics.babraham.ac.uk/projects/fastq_screen>), and visualized with MultiQC (v0.8). Adaptor sequences, first 9 bp and low quality ends were removed with Trim Galore (v0.4.2, <http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/>; parameters:--phred33--length 36--clip_R1 9). Trimmed reads were aligned to a combined human (GRCh38) and mouse (GRCm38) genome reference using STAR (v 2.5.2a; parameters:--runMode alignReads--genomeLoad LoadAndExit) ([@bib13]) in order to exclude potential contaminations in the xenograft models. Quality control of resulting BAM files was performed with RSeQC (v2.6.4) ([@bib65]). The gene-level expression counts were computed with featureCounts (v1.5.1) ([@bib41]) using the gene annotations obtained from RefSeq (release 70). Samples with less than 800 features detected (threshold ≥ 1 mapped read) or showing more than 5% of contamination from the other species were removed from further analysis. To normalize gene counts for cell-specific biases in single-cell data, we used size factors computed using the deconvolution strategy developed by [@bib44] and implemented in the scran (v1.6.5) package available on R/bioconductor ([@bib24]). RNA-seq data processing, quality control and visualization was performed with the help the R/bioconductor package scater (v1.6.0). Gene regulatory network analysis to identify putative transcription factor dependence was generated using iRegulon ([@bib25]).

### Weighted Gene Co-Expression Analysis {#sec4.5.6}

A signed co-expression network was constructed using Weighted Gene Correlation Network Analysis (WGCNA) R package (v1.61) based on 335 genes from a well-established stem cell gene signature ([@bib67]). The analysis was performed only on those genes with a Fragments Per Kilobase Million (FPKM) value ≥ 1 in at least three samples (n = 301). A soft-threshold power used to create the pairwise distance matrix was defined as the lowest power for which the scale-free topology fit index reaches 0.80. Dynamic tree cut was (dynamicTreeCut R-package, v1.63) used to defined modules of genes. The association between the modules eigenvalues and CTC groups was evaluated using Wilcoxon rank sum. Gene ontology network analysis was generated using ClueGO ([@bib63]).

Data and Software Availability {#sec4.6}
------------------------------

### Data Analysis and Data Availability {#sec4.6.1}

Data analysis, statistical test and visualization were conducted in R (version 3.4.0; R Foundation for Statistical Computing, Vienna, Austria) and GraphPad Prism (v7.0). RNA sequencing data have been deposited to Gene Expression Omnibus ([@bib14]). The accession number for the RNA sequencing data reported in this paper is GEO: [GSE111065](ncbi-geo:GSE111065){#intref0075}. Bisulfite sequencing data have been deposited to the European Nucleotide Archive. The accession number for the bisulfite sequencing data reported in this paper is ENA: PRJEB25101.

Supplemental Information {#app2}
========================

Table S1. Clinical Features of Breast Cancer Patients Prior to CTC Collection, and Number and Features of Single CTCs and CTC Clusters Isolated from Patients and Xenografts, Related to Figure 1Table S2. Whole-Genome Bisulfite Sequencing Sample Statistics, Related to Figures 1 and 6Table S3. RNA Sequencing Sample Statistics, Related to Figures 3 and 6Table S4. WGCNA of Genes (≥3 TPM Counts) Significantly Enriched (p \< 0.05) in Patient- and Xenograft-Derived CTC Clusters, Xenograft-Derived Single CTCs, as well as Stemness-Related Genes in Patient- and Xenograft-Derived CTCs, Related to Figures 3 and 4Table S5. Expression Modules Enriched in Patient-Derived CTC-Clusters (Blue, Gray), Single CTCs (Turquoise, Brown) and Xenograft-Derived CTC Clusters (Green, Yellow), Related to Figure 4Table S6. FDA-Approved Compounds Used in the Screen, Related to Figure 5Table S7. Stemness-Related Genes Up- or Downregulated in BR16 and BRx50 as a Result of a 17-Day Treatment with 20 nM of Digitoxin or Ouabain, Related to Figure 6

We thank all patients that participated in our study, all involved clinicians and study nurses, and all members of the Aceto lab for feedback and discussions. We thank Dr. Joan Massagué (MSKCC, New York); Dr. Daniel Haber and Dr. Shyamala Maheswaran (MGH and Harvard Medical School, Boston, MA) for donating cell lines; and Dr. Gerhard Christofori (University of Basel, Switzerland) for helpful comments on the manuscript. We thank Katja Eschbach and Elodie Burcklen from the Genomics Facility Basel (ETH Zürich) for sequencing support; Falco Kilchmann and Eva Maria Riegler from NEXUS PHT (ETH Zurich) for help with the drug screening; and Dr. Robert Ivanek (University of Basel) for support with data analysis. Calculations were performed at sciCORE scientific computing center at the University of Basel. Research in the Aceto lab is supported by the European Research Council (678834), the Swiss National Science Foundation (PP0P3_163938), the Swiss Cancer League (KFS-3811-02-2016 and KLS-4222-08-2017), the Basel Cancer League (KLbB-4173-03-2017), the two Cantons of Basel through the ETH Zürich (PMB-01-16), and the University of Basel.

Author Contributions {#sec5}
====================

S.G. and N.A. conceived the study, generated and interpreted the results, and wrote the manuscript. F.C.-G. conceived and performed the bioinformatics analyses. B.M.S. generated RNA sequencing (RNA-seq) libraries from CTCs. M.V., J.L., C.K. C.R., V.H.-S., and W.P.W. provided clinical input, consented breast cancer patients, and provided blood samples for CTC isolation. R.S. and I.K. performed cell culture experiments and contributed to the HTS compound screen. M.C.S. performed the CRISPR-related experiments. C.B. performed the WGBS and the RNA-seq library sequencing. C.U.S. provided the FDA-approved library for the drug screen. All authors approved the manuscript.

Declaration of Interests {#sec6}
========================

N.A., S.G., and I.K. are listed as inventors in patent applications related to cancer treatment through CTC clusters inhibition. All other authors declare no competing interests.

Supplemental Information includes seven figures and seven tables and can be found with this article online at [https://doi.org/10.1016/j.cell.2018.11.046](10.1016/j.cell.2018.11.046){#intref0080}.

[^1]: Lead Contact
