Expanded genomic profiling of circulating tumor cells in metastatic breast cancer patients to assess biomarker status and biology over time (CALGB 40502 and CALGB 40503, Alliance)

Patient samples and CTC enumeration

Clinical samples were collected between May 2010 to May 2012 from MBC patients enrolled in local (University of California San Francisco, UCSF) and nationwide clinical trials Cancer and Leukemia Group B (CALGB) 40502 (NCT00785291) and 40503 (NCT00601900), and Translational Breast Cancer Research Consortium (TBCRC) 009 (NCT00483223) (Supplementary Table 1). CALGB is now part of the Alliance for Clinical Trials in Oncology. All patients gave informed consent under a protocol approved by institutional review boards in respective participating centers. Blood was drawn into CellSearch^™ CellSave preservative tubes (Veridex) for CTC enumeration (Figure 1a). Additional volumes of blood were drawn into tubes containing EDTA for CTC isolation. Blood samples from other institutions were shipped overnight to the Park Laboratory at UCSF and processed immediately for CTC enumeration. Previous studies have shown that CTC numbers obtained from CellSearch^™ and IE/FACS (described below) were highly concordant (2, 10). To determine whether a blood sample contained CTCs for isolation, we first screened blood samples using the CellSearch^™ assay. Blood samples with ≥8 CTCs/7.5mL (~ ≥1 CTC/mL) were subjected to IE/FACS. CTC isolation was performed between 24–36 hours after blood draw. Both CTC enumeration (CellSearch) and isolation (IE/FACS) were performed in the Park Laboratory at UCSF. Detailed descriptions of the methods are provided in the Supplementary Information.

Cell isolation by IE/FACS

CTCs were isolated via IE/FACS as previously described (4–6) (Figure 1b). Briefly, magnetic beads coated with EPCAM monoclonal antibody were used to enrich for tumor cells. The tumor-enriched samples were then subjected to FACS analysis using differentially labeled monoclonal antibodies to distinguish tumor cells (nucleated, EPCAM+/CD45−) from leukocytes (nucleated EPCAM−/CD45+) during cell sorting. IE/FACS allowed for the isolation of single or small pools of cells, which could then be subjected to downstream molecular analyses.

Panel of 64 genes for transcript analysis

We used a panel that has been previously validated for gene expression analysis of rare tumor cells (11). The list includes EPCAM and PTPRC (encodes CD45), which are markers for epithelial cells and leukocytes, respectively. Also included are clinically relevant cancer genes (e.g., ESR1 and ERBB2), stem cell and epithelial-mesenchymal transition (EMT) markers, and candidate reference genes for data normalization (Fig. 1C). The complete list can be found in Supplementary Table 2.

Taqman Low-Density Array Quantitative PCR (aQPCR) analysis

For multiplex aQPCR analysis, we used a custom microfluidic card (384-well format) containing two sets of 64 Taqman probes printed in triplicate. Details of the cell lysis, reverse transcription, specific transcript (cDNA) preamplification, and aQPCR analysis are described in the Supplementary Methods. Results of performance evaluation of the aQPCR assay are presented in the Supplementary Results and Supplementary Figure 1a–e.

Reference genes

To select the reference genes for normalization, we used the geNorm algorithm within RealTime StatMiner^® (see below) to calculate the gene stability measure (M) for candidate genes: ACTB, GAPDH, GUSB, RPLP0, TFRC, and RPS18. Genes ACTB and RPS18 were chosen because they showed the lowest M values indicating least variable expression across 105 CTC and 76 leukocyte samples (Supplementary Figure 1a).

Quality controls

Prior to aQPCR analysis, preamplified cDNA samples were screened for expression levels of ACTB, GAPDH, and RPS18 using conventional QPCR. Results of initial experiments revealed that samples with low transcript levels of RPS18 (Cycle threshold, Ct ≥26) (12) and/or ACTB and GAPDH (Ct >36) resulted in a failed aQPCR analysis, indicating insufficient quantity and/or poor quality RNA. Samples with <20% detection rate (i.e., detection of ≤12 of the 64 genes by aQPCR analysis) and those with undetectable expression (missing values) of references genes, RPS18 and ACTB, were excluded from the analysis.

Single cell gene expression analysis

Single cells were sorted by IE/FACS and subjected to expression profiling. Transcript levels for a subset of 32 genes were measured via the Fluidigm^™ method following manufacturer’s instructions (Supplementary Table 2). Details of the protocol are described in the Supplementary Methods.

Array comparative genomic hybridization analysis (aCGH)

CTCs collected from a single time point were analyzed by aCGH in 49 of the 102 patients in the study. Genome-wide copy number aberrations were assessed using bacterial artificial chromosome aCGH, as previously described in detail (5, 13). Microarray data was subjected to circular binary segmentation using the DNAcopy package in Bioconductor (14, 15). Details of the computational methodologies for copy number analysis including the use of the Nexus 6.0 software (Biodiscovery) are described in the Supplementary Methods.

Gene expression analysis

We used the RealTime StatMiner^® (version 4.2) to analyze gene expression data as previously described (11). The mean Cts for ACTB and RPS18 were used to calculate the ΔCts. Cluster analysis was performed using unsupervised ward linkage hierarchical clustering methods with Pearson correlation distance as a similarity measure. Multidimensional scaling analysis was also performed to visualize similarities in gene expression profiles among samples. Differentially expressed genes between two groups were elucidated using a parametric analysis (Limma) for unpaired samples, and a paired T-test for paired samples. To adjust for multiple comparisons, the Benjamini-Hochberg (BH) correction method was used (16), and an adjusted p value of <0.001 was considered statistically significant. Relative quantification (RQ) was reported in the logarithmic scale (log₁₀RQ=log₁₀ 2^−ΔΔCt). A log₁₀RQ=0 indicated no differential expression, log₁₀RQ=1 or −1 indicated that a gene is expressed 10 times or 1/10 as much in the test sample relative to the calibrator sample, respectively. Percent detection for each of the 64 genes in CTC and leukocyte samples was compared using a two-tailed Z-test. A p value of <0.05 was considered statistically significant.

Serial expression analysis

To determine whether CTC samples from the same patients exhibit a unique gene expression signature, we calculated the mean Pearson and Fisher-transformed Pearson correlation (17, 18) as well as the mean Euclidean distance between samples from the same patient and for all other patients. Patient IDs were then randomly permuted 10,000 times and the mean distance between samples per ‘simulated patient’ was calculated. A p-value was calculated using the Fisher transformation (17) for actual and simulated random data to test whether serial CTC expression profiles from the same patients were more similar to one another compared to those from other patients.

Assessment of ESR1/ERBB2/MKI67 status and intrinsic subtype analysis

We determined ESR1 status in CTCs using gene expression data. The bimodal distribution of ESR1 expression allowed for dichotomization into positive and negative groups using the local minimum between modes as a cutoff point. For ERBB2 status, CTCs were classified as ERBB2-positive based on aCGH defined copy number gains or amplification of the ERBB2 locus. For CTC samples with no aCGH data (n=50), ERBB2 status was dichotomized based on a simple aQPCR gene expression-based classifier; this was derived by receiver-operating characteristics curve analysis (Youden index) of a subset of samples with both gene expression and aCGH data (n=101) (Supplementary Figure 1f). For proliferation status, we used MKI67 mRNA expression as a surrogate marker for Ki67 protein expression. We defined the top tertile (threshold=66% quantile) as ‘high’ and the remainder as ‘low’.

Intrinsic subtype was determined using the 16 PAM50 classifier genes present in the assay. We used a previously validated program and algorithm described by Parker and colleagues to make subtype assignments (19). Subtype calls were filtered using confidence measures produced by the program, only samples with >80% confidence were reported (i.e., 69/105 (66%) of patient CTC samples). To validate our approach, we compared the intrinsic subtype calls using all of the PAM50 genes vs. the subset of 16 using the dataset available at https://genome.unc.edu/pubsup/breastGEO/PAM50.zip, and observed 80% concordance (20).

Analysis and visualization of the data was performed using packages in R (21).

Survival analysis

We examined association between with CTC profiles and patient outcomes. The endpoints included progression-free survival (PFS), defined as time from study entry to documented progression or death from any cause, and overall survival (OS), defined as time from study entry to death from any cause. PFS and OS were estimated using the Kaplan-Meier method (22), and the log-rank test (23) was used to compare survival between groups, with p<0.05 being considered statistically significant. The survival data used to assess outcome were collected by the Alliance Statistics and Data Center. Statistical analyses were performed on all data available as of July 18, 2016.

Development and performance testing of CTC assays

DNA copy number profiling of CTCs was performed using our previously validated approach of whole genome amplification of CTC genomic DNA followed by aCGH analysis (5, 11).

For gene expression analysis, we developed a protocol involving microfluidic-based multiplex QPCR array (aQPCR) to analyze the transcript levels of 64 cancer-related genes. The method involves reverse transcription of total RNA, preamplification of specific cDNA transcripts and aQPCR analysis. A detailed discussion of the assay optimization and performance analysis is presented in the Supplementary Materials and results are displayed in Supplementary Figure 1b–e.

Proof-of-concept experiments using spike-in models

The performance of IE/FACS and aQPCR assay was further evaluated using spiked cancer cell lines (n=6) into healthy blood or cell admixtures (Supplementary Figure 2). Gene expression analysis of captured cells revealed high purity and specificity of IE/FACS isolation, and the robustness of the preamplification/aQPCR protocol (see Supplementary Materials for details).

Analytical validation of the CTC gene expression assay in clinical samples

We applied our optimized protocols to characterize CTCs from a prospective series of MBC patients who had provided informed consent in different multicenter clinical trials. The workflow of the study is outlined in Figure 1a. Starting from 244 blood samples from 162 consecutive patients, we isolated ~10–20 CTCs by IE/FACS and performed gene expression analysis (Figure 1b–c, Supplementary Table 1). Of the 244 samples, 93 samples were excluded due to low CTC yield (n=46), poor RNA quality (n=37), or failed aQPCR experiment (n=10) (Supplementary Figure 3). The remaining 151 (61%) samples from 105 unique patients yielded high quality expression data and were used for subsequent analysis. The dataset also included additional 46 CTC expression profiles collected at later time points from 28 patients (see below).

Leukocytes were isolated from the same blood samples, and profiled in parallel as non-tumor controls. Matched leukocyte expression data was available for 76 of the 105 patients, 29 patients lacked data because the samples had poor quality RNA.

Cluster analysis showed that CTCs grouped among themselves and away from normal blood (Figure 1d). Multidimensional scaling analysis to visualize expression data from all CTC (n=151) and leukocyte samples further confirmed that CTC expression profiles were clearly distinguishable from those from those of leukocytes (Figure 1e)

Transcripts coding for EPCAM, CD24, KRT19, ERBB2, CCND1, CAPG, BAG1, as well as reference genes, including ACTB, RPS18, RPLP0, GAPDH, and GUSB, were frequently detected (>90%) in CTCs (Supplementary Figure 4a).

Expression of EMT, stem cell, and other cancer-related genes in CTCs

To determine differentially expressed genes between CTCs and leukocytes in 76 paired samples (BH corrected p<0.001), we compared the relative expression (ΔCt) of genes between the two groups with leukocytes as the calibrator. Consistent with our epithelial cell isolation approach strategy, CTCs exhibited up-regulation of EPCAM and down-regulation of PTPRC/CD45 (Figure 1f). Other breast cancer-related genes including CCND1, KRT7, MUC1, and TFF3 were up-regulated while CD68, a monocyte/macrophage marker, was down-regulated. Expression of EMT (SNAI1: 51%, TWIST1: 22%, CAV: 24%, and VIM: 78%) and stem cell markers (CD44: 88%, CD24: 100%, and ALDH1A1: 22%) were observed in a subset of CTCs (Supplementary Figure 4b). However, when compared to leukocytes, all four EMT markers were down-regulated in CTCs, as were stem cell markers CD44 and ALDH1A1; while CD24 was up-regulated (Supplementary Figure 4c). Taken together, our results with isolated EPCAM-positive CTCs do not indicate epithelial/mesenchymal or stem cell-like phenotypes.

Receptor and intrinsic subtype analysis

We classified CTC samples from 105 patients into groups according to receptor and intrinsic subtypes. ESR1 displayed a bimodal distribution allowing dichotomization into positive and negative groups (Figure 2a). MKI67 mRNA expression, which also showed a bimodal distribution, was used as a surrogate proliferation marker for Ki67 expression (Figure 2b). For ERBB2 status, we used both copy number and expression data (see Methods).

Approximately 70% of the CTC samples were considered ESR1-positive, of which 48% were ESR1+ERBB2− and 22% ESR1+ERBB2+ (Figure 2a). ERBB2-positive (and either ESR1-positive or negative) and ESR1−ERBB2− CTCs accounted for 27% and 26% of the samples, respectively. 65% of the CTC samples were considered to have low proliferative (MKI67) status. Of the 105 samples, 69 were assigned intrinsic subtype calls using the 16 PAM50 genes present in the assay (see Methods). These samples were subdivided into 30% luminal A, 6% luminal B, 13% HER2-enriched, 33% basal, and 12% normal-like subtypes. As expected, we found a significant association between intrinsic and receptor subtype calls (Supplementary Figure 5a, chi-squared test, p=0.00055).

Clustering analysis using the 16 PAM50 genes revealed three major clusters (Figure 2b). The leftmost cluster (Cluster 1) included CTC samples comprising mostly basal and/or ESR1−ERBB2− profiles with high proliferative status. This cluster contained CTCs with matched primary tumors that were also ER-negative. Luminal B CTCs were found only in Cluster 2, while luminal A CTCs were found in both Clusters 2 and 3, and so were ESR1+ERBB2− and ESR1+HER+ receptor subtypes. As expected, similar results were obtained when all 64 genes were used for clustering analysis (Supplementary Figure 5b).

Monitoring CTC-based biomarker status over time

To assess the feasibility of serial CTC expression analysis over time, we isolated and profiled CTCs in 74 serial blood samples from 28 patients. Blood was collected from each patient 2 to 5 times at intervals ranging from 17 to 307 days from the first sample collection (T1). Visual inspection of the expression profiles of CTC samples from individual patients showed fluctuations in expression at the individual gene level. However, intrinsic and clinical biomarker based phenotype assignments were generally consistent across time points (Figure 3a, Supplementary Figure 6a). 100% (28/28) of patients who were assessed for CTCs at multiple time points showed consistent ERBB2 status (by aCGH/aQPCR) at all time points sampled. ESR1 and proliferation (MKI67) status (by aQPCR) were less consistent (2-sample proportion test: p=0.0036 and p>0.001, respectively), with 32% (9/28) of patients showing a change of ESR1 status over time and 54% (15/28) exhibiting variable proliferation status.

To quantify changes in gene expression in serial CTC samples, we calculated the correlation between the first two samples (T1 vs. T2) from each patient (n=28), while noting the time interval between the two blood draws. We observed that CTCs collected and analyzed further apart in time—i.e., greater number of days elapsed between T1 and T2—showed lower Pearson correlation than those with shorter intervals between sample collections (r=−0.434, p=0.021, Figure 3b, Supplementary Figure 6b). To extend this analysis, we examined the correlation of T1 samples with all samples from individual patients, including those collected at later time points. We observed that, in general, the median Pearson correlation of patients’ samples decreased across time points, although this pattern was not statistically significant (Figure 3c, Supplementary Figure 6c). Using the information above, we estimated an average Pearson correlation change of −0.31/year.

Visual examination of the expression profiles showed similarities in gene expression in serial CTC samples from individual patients. To test whether serial CTCs from a particular patient (ranging from 2 to 5 samples) were more similar to one another than to CTCs from other patients, we compared the correlation between metachronous samples from individual patients vs. correlations between randomly paired samples from a large simulated dataset (generated by permuting sample labels 10,000 times). This comparative analysis revealed that correlation among samples from the patient data (Fisher transformed mean ρ=0.66) were significantly higher compared to those from the simulated dataset (Fisher transformed mean ρ=0.80), indicating that CTC samples from individual patients were more related to each other than to randomly selected samples (p<0.0001; Figure 3d, Supplementary Figure 6e). Moreover, visualization by multidimensional scaling analysis in four patients with the highest number of serial samples showed that CTCs from the same patient over time are more similar to one another than to samples from other patients (Figure 3e). Taken together, these observations indicate the reproducibility of our assay, and suggest that CTC expression profiles can be obtained from individual patients.

Single cell gene expression analysis reveals heterogeneity among individual CTCs

To assess tumor heterogeneity at the single cell level, we performed expression profiling on individual CTCs by multiplex aQPCR analysis using a microfluidic-based dynamic array (Fluidigm^™). Pooled and single CTCs were isolated using IE/FACS from two MBC patients: 4042 (Figure 4a) and 4043 (Figure 4b). Expression analysis revealed higher levels of heterogeneity observed among single cells compared to pooled cells (Figure 4c). Despite the heterogeneity, some genes (e.g., ESR1) did exhibit more consistent expression across all samples. Multidimensional scaling (Figure 4d) and clustering analysis (Figure 4e) showed that CTCs from each patient displayed a unique gene expression signature, consistent with results from serial analysis of CTC samples. In addition, the expression signature seemed to be conserved down to individual cells, as single cell samples from each patient clustered together (Figure 4d).

Copy number analysis of CTCs reveals aberrations frequently seen in primary breast cancers

Duplicate pools of CTCs were isolated from the same tumor-enriched blood samples in 49 of the 105 patients for aCGH profiling. Genome-wide copy number analysis detected genomic aberrations in all of our CTC samples (Figure 5a; Supplementary Figure 7a). Frequent aberrations included gains in 1q and 8q, as well as losses in 8p and 16q (Supplementary Table 3). Gene amplifications were also detected for oncogenes, such as CCND1 (31%), ERBB2 (12%), and MYC (29%). Similar recurrent aberrations were observed when compared to a publicly-available dataset for primary breast tumors, consistent with our previous findings (Supplementary Figure 7b) (24).

Cluster analysis of copy number data revealed three distinct groups (Figure 5b). Examples of CTC copy number profiles for each cluster are shown in Supplementary Figure 7c. Cluster 1 contained CTCs that exhibited significantly less genomic aberrations compared to CTCs in Clusters 2 and 3 (Supplementary Figure 7d). Cluster 2 included CTCs with 8q gain, while Cluster 3 included those with 1q gain/11q loss. Interestingly, CTCs with ESR1−ERBB2−/basal phenotype were mostly observed in Cluster 1 containing CTC samples that display low genomic instability. Associations between genomic aberrations detected in CTCs vs. patient outcomes were not analyzed due to the limited sample size.

ERBB2 copy number and expression in CTCs

Our previous studies have shown that CTC aCGH analysis can accurately detect copy number changes involving ERBB2 (5, 11). We therefore performed aCGH analysis of the CTCs in the present study, and observed ERBB2-positive and ERBB2-negative CTCs (Figure 5c).

Other studies have shown that ERBB2 mRNA expression and copy number in primary breast cancer are highly concordant (25, 26). To evaluate concordance in CTCs, we compared ERBB2 expression and copy number from 49 patients in this study and an additional 88 patients from our previous work (5). Comparative analysis showed that the mean expression levels of ERBB2 in ERBB2-positive CTCs (by copy number) was numerically higher compare to ERBB2-negative CTCs; however, the agreement between copy number and expression data was not strong (Figure 5D). Also, broad ranges of ERBB2 expression levels was observed in both ERBB2-positive and ERBB2-negative CTCs (by copy number), but more so in the latter group.

Biomarker status in CTCs and matched primary tumors

To investigate whether biomarkers relevant to breast cancer can change during disease progression, we compared the ESR1 (by aQPCR) and ERBB2 (by aCGH) status in CTCs with clinical ER and HER2 status of matched primary tumor. Comparative analyses between CTCs and corresponding metastases were not performed because metastatic tissue was not available. Also, we did not analyze gene expression in matched primary tumors.

Of the 105 patients in this study, 70% harbored ESR1-positive CTCs. ER status for 102 corresponding primary tumors was available for comparison, of these 85% were ER-positive. ER status was in agreement in 73% (74 of 102) matched samples (Figure 5e). Conversely, 25% (22 of 87) of patients with ER-positive primary tumors showed ESR1-negative CTCs, while 40% (6 of 15) of patients with ER-negative primary tumors had ESR1-positive CTCs.

Assessment of ERBB2 status using available aCGH data revealed that CTCs from 26% (35 of 137) of the patients showed copy number gains, and therefore were considered ERBB2-positive. Comparison of ERBB2/HER2 status between 130 matched CTCs and primary tumors revealed a concordance of 77% (100 of 130) (Figure 5f). Of note, ERBB2-positive CTCs were detected in 24% (30 of 127) of HER2-negative primary tumors.

CTC profiles vs. patient outcome

We examined the correlation between CTC biomarker status and patient outcome. Follow-up period was approximately 40 months. There were no significant differences in OS between patients with ESR1-negative CTCs vs. those with ESR1-positive CTCs (Supplementary Table 4, Supplementary Figure 8a), and patients with ERBB2-negative CTCs vs. those with ERBB2-positive CTCs (Supplementary Figure 8b). We did observe that patients harboring CTCs with a more aggressive subtype (i.e., ESR1-negative or ERBB2-positive) had numerically shorter survival compared to those with less aggressive phenotype.

Survival analysis based on receptor subtype (combined ESR1 and ERBB2 status) revealed that patients with ESR1−ERBB2+ CTCs have significantly shorter PFS (log-rank p=0.015) compared to other subtypes (Supplementary Figure 9a). No significant differences were observed for OS (log-rank p=0.68). Patients with basal-like CTCs showed a trend towards shorter PFS (log-rank p=0.1053) and OS (log-rank p=0.1322), which was not statistically significant (Supplementary Figure 9b). Analysis based on proliferation status revealed that patients with CTCs showing high expression of MKI67 (encodes the proliferative marker Ki67) have significantly shorter PFS (log-rank p=0.0011, Figure 6a) and OS (log-rank p=0.01, Figure 6b) compared to patients with CTCs that express the gene at low levels.

ESR1 status in CTCs and ER status in primary tumors vs. patient outcome

We also investigated whether combining biomarker status in primary tumors and CTCs was correlated with patient outcome. Stratification based on ER status revealed that patients whose CTCs and primary tumor were both ESR1−/ER− had significantly shorter PFS (log-rank p=0.01, Figure 6c) and OS (log-rank p=0.0035, Figure 6d) compared to patients who were ESR1+/ER+ in both compartments. In patients where ESR1/ER status in primary tumor and CTCs changed from positive to negative or vice versa, the PFS and OS were similar to survival probabilities based on the original ER status of the primary tumor. Survival analysis based on ERBB2/HER2 status was not performed due to the limited number of evaluable HER2-positive patients in the study.