To better understand heterogeneity at the single-cell level, we orthotopically implanted four patient-derived xenografts (PDX) [8, 9] and a cell line (see Supplementary Table 1) known for their lung metastatic potential into NOD-SCID-Il2rgnull (NSG) mice (Fig. 1a). Tumours were resected from the mammary fat pad and the animals were monitored for metastasis. Once the animals showed signs of distress (i.e., weight loss, difficulty to breath), lungs presenting metastatic lesions were collected and processed for single-cell transcriptional profiling. To exclude murine cells from the downstream analysis, human cancer cells from the tumours and matched lung metastases were purified via FACS gating GFP-positive MDA-MB-231 cells (Fig. 1b) or CD298-positive cells for PDX models (Fig. 1c). Single cells were isolated using a microfluidic device (Fluidigm C1) ahead of library preparation and sequencing. This workflow yielded a total of 1,523 single cells (Supplementary Fig. 1a) after RNA-Seq and quality control.

Interpatient heterogeneity is dominant over intrapatient heterogeneity. a Overview of the experimental setting and models used. Human breast cancer models were implanted in the mammary fat pad of NSG mice. Tumours and lung macrometastases were harvested and mechanically and enzymatically dissociated. Human cells were purified by FACS using GFP or CD298 staining. Single cells were isolated with the Fluidigm C1 microfluidic platform and then sequenced. b Representative FACS strategy for the isolation of MDA-MB-231 GFP positive. d tSNE plot showing initial clustering of all the sequenced cells, according to models of origin. e tSNE plot showing initial clustering of all the sequenced cells, according to the site of origin

Initial clustering of the quality-controlled data revealed that the cells formed groups (clusters) according to the donor models (Fig. 1d). Within each cluster, cells did not clearly separate based on their origin, tumour, or metastasis (Fig. 1e). These observations highlight the importance of interpatient over intrapatient heterogeneity.

We then asked whether cells gather within each donor-specific cluster by known biological or technical features. We projected the percentage of detected genes (PDG) in each single-cell library (Fig. 2a) onto t-distributed stochastic neighbor embedding (t-SNE) [10] and found that the PDG influences the clustering of the cells within each given model. This variable potentially represents both biological and technical effects. Next, we assessed whether the cell cycle stage influenced the analysis of the data, as this biological variable has a broad impact on gene expression [11]. We developed a method to infer cell cycle stages in single cells (in R package gripgh) and each cell was labelled with one of the four labels (G1, G1/S, S/G2, G2/M). We found that the cell cycle stage does influence the clustering of the cells within each model (Fig. 2b).

Cell cycle and percentage of detected genes delineate cell clustering. Single cells within PDX models cluster by a, library complexity (PDG) and b, cell cycle similarity. c Clustering performed via t-SNE produced 16 cell clusters. d Cluster composition according to cell cycle (left) and library complexity (right) e An in-silico cell cycle scoring prediction to arrange single cells on a continuous cell cycle spectrum in addition to distinct cell cycle stages

We then applied graph-based clustering to the different models and obtained 16 subclusters (Fig. 2c). These subclusters were mostly composed of cells in a similar cell cycle stage (Fig. 2d, left bar graph) and with similar PDGs (Fig. 2d, left bar graph). The influences of both the cell cycle and the PDG were also observed when the models were analysed individually (Supplementary Fig. 2a, top bar graphs for each model). Marked interpatient heterogeneity also led to the clustering of the cells according to models (Supplementary Fig. 2b). Altogether the data suggest that, in contrast to the site of origin (i.e., primary tumour or metastasis), the PDG and the cell cycle both influence the clustering of cells (Fig. 1e).

For the cell cycle stage prediction we considered that the cell cycle is not composed of discrete stages but is more a continuum of states. Gene expression is gradually modulated as a cell progresses in the cycle. Using this riche single-cell RNA sequencing data and genesets whose expression varies in different cell cycle stages, we created a circular trajectory and placed the cells in a cycle (Fig. 2e). This precise allocation along the cell cycle continuum allowed precise cell cycle staging.

To observe underlying biological processes involved in different tumour cells and to (beyond donor effect, PDG and cell cycle stage) and group biologically similar cells, we needed to remove the confounding factors of donor effect, PDG and cell cycle stage. To remove biases attributed to these factors, we used a generalized linear model (GLM) [12], correcting gene expression according to the position of the cell on the cell cycle spectrum and the complexity of the RNA-Seq library that it yielded (Fig. 3a).

Removal of cell cycle variation and percentage of detected genes reveals three major biological clusters. a Generalized linear model approach used to remove biases (cell cycle, and library complexity). b Post-correction the new cell clusters have a more balanced distribution of cells from different cell cycles, library complexity (PDG), and cell source (PDX model). c Gene Set Enrichment Analyses (GSEA) of individual clusters reveal common and different Hallmark genesets enriched among clusters. Top 25 of these genesets according to the absolute values of the normalized enrichment score (NES) are shown and were used to identify the superclusters

Initially we performed gene set enrichment analysis (GSEA) with the Hallmarks gene set [13] (Supplementary Fig. 3a) on clusters defined for each model on the corrected data (Supplementary Fig. 2a, bottom bar graphs for each model, named Ax, Bx, Fx, Dx, Ex). We observed that cell clusters from different models show enrichment in a common set of biological processes. This suggests that each model contains cells in closely related biological states that could now be visualized after bias correction.

We then analysed the cells irrespective of their models or sites of origin. Using the corrected gene expression data, the cells formed 14 new clusters (Cx). The compositions of these groups were then analysed according to cell cycle stage (Fig. 3b, left), PDG (Fig. 3b, centre), and model (Fig. 3b, right). Correcting for these variables led to unbiased clustering of the cells, with clusters composed of cells from various cell cycle stages, library complexities, and models. We also plotted the composition of each cluster in terms of site of origin before and after correction (Supplementary Fig. 3b). Corrected clusters exhibited a more balanced composition, with roughly equal proportions of cells originating from the tumour and lung metastases. The data suggest that populations of cells clustering together due to their biological similarities can be found in both the tumour and the metastatic sites, without specificity to one or the other.

To further investigate the different biological states suggested by the data in Supplementary Fig. 3a, we performed GSEA on the 14 clusters formed by all the cells, regardless of their model of origin. The 14 clusters formed three “super” biological clusters (Fig. 3c). Supercluster A (C8 to C6, Fig. 3c left) is characterized by low enrichment for most of the Hallmark geneset pathways. Supercluster B (C13 to C2, Fig. 3c centre) displays the most heterogeneous regulatory landscape, with highly, moderately, and minimally enriched processes. This supercluster is defined by highly enriched Hypoxia and TNFα signalling via NFκB. The moderately enriched pathways include relevant processes such as EMT, TGFβ signalling, Interferon Gamma response, or the P53 Pathway. Finally, the least enriched genesets include MYC signalling, G2M checkpoint, and E2F targets. Interestingly these processes are most enriched in supercluster C (C7 to C4, Fig. 3c right), which also displays marked enrichment of genesets pertaining to oxidative phosphorylation, mTOR signalling, fatty acid metabolism, and DNA repair. Genesets enriched in supercluster B suggest a phenotype related to EMT, while cells in supercluster C appear to be proliferating while still partially enriched for EMT-related pathways.

We then plotted the repartition of the cells according to supercluster allocation (Supplementary Fig. 4a), noticing that cells from superclusters B and C were the most distant, with cells from supercluster A in between. We also assessed repartition according to the site of origin (Supplementary Fig. 4b) and once again found a relative balance between tumour origin and lung metastases origin of the cells forming the superclusters.

Next, we selected EMT and Proliferation markers significantly altered (FDR < 0.05) in pairwise comparisons between the superclusters (Fig. 4a). Proliferation markers Ki67, MCM3, and PCNA [14] confirmed that supercluster C is the most proliferative, followed by cluster A, while supercluster B expresses these markers the least. EMT markers indicated that this process was taking place at varying levels across the different superclusters. Supercluster A displayed the least engagement in the transition according to its low expression of several EMT markers (ZEB1, SOX9, SNAI1, FN1, TGFBR1). Superclusters B and C showed increased expression of these markers but at varying levels. Such heterogeneity suggests that these superclusters may undergo EMT but could be at different stages of the process. Such partial EMT has previously been described [15, 16] and may reflect the balance between proliferative potential and migratory capability, both properties being typical of different stages of the metastatic cascade.

Comparison of major biological superclusters reveals partial EMT state regulators. a EMT and Proliferation markers significantly altered (FDR < 0.05) in supercluster pairwise comparisons. b, c and d, Right: Transcription Factors Enrichment Analysis of the top 250 up- and downregulated differentially expressed genes of each supercluster. Left: Volcano plots highlighting the cell cycle gene set (blue, Bioplanet 2019), the EMT gene set (orange, mSigDB) and hits of interest (red) in corresponding supercluster pairwise comparisons

To investigate partial EMT states of the superclusters, we selected the top 250 up- and downregulated genes of each cluster and performed GSEA as well as Transcription Factor Enrichment Analysis (TFEA) with the EnrichR [17] platform (Fig. 4a, b, c bar plots).

The “TF Perturbations followed by expression” geneset, which was generated by the curation of experiments altering TFs before measuring gene expression, revealed that different TFs govern the top differentially expressed genes (DEGs) across the superclusters. These uncommon TFs paint a complex picture of the EMT states existing in superclusters A, B, and C. HEY2, OVOL2, TFAP2C, and TP63 are TFs described as regulators of partial EMT states in mouse models [16]. ZNF750 was recently described as an EMT repressor in breast cancer [18], an activity shared by OVOL2 [19], TP63 [16], and FOXO1 [20]. HEY2, TFAP2C  [21, 22], SOX4, and SOX9 [23] are known to promote EMT. U2AF1 is a splicing factor fine-tuning translation with reported effects in development and EMT.

These TFs control the EMT and proliferation state of the superclusters shown in Fig. 4a. Supercluster A, shown to be mildly proliferative, is controlled by TFs evocating differentiated slowly proliferating cells. E2F4 is known to be abundant in differentiated cells and to repress proliferative genes [24]. TP63 and TFAP2C have been described as controlling early hybrid EMT states, with cells close to an epithelial state and more prone to proliferate than mesenchymal cells. Additionally, KDM5B has been reported to characterize a slow-cycling cell subpopulation in melanoma [25], which fits the traits of supercluster A. Supercluster B, the least proliferative, is also the one with the strongest expression of canonical EMT markers/TFs such as FN1, SNA1, MYC, and SOX9. When compared to superclusters A and C, DEGs from supercluster B appear to be under the control of HEY2, a member of the basic helix-loop-helix (bHLH) transcription factor family. bHLH have been described as late hybrid EMT regulators [16] responsible for mesenchymal phenotypes that favour migration and have little proliferation potential. Supercluster C, the most proliferative, is controlled by ZNF750, TFAP2C, and OVOL2. These TFs have been described as either EMT repressors or regulators of early hybrid EMT stages, corresponding to epithelial-like phenotypes permitting proliferation. It should however be noted that the EMT markers are high in this supercluster, indicating that the process is ongoing yet likely oriented towards the proliferation of cells rather than their migration.

Next, we analysed individual hits from the pairwise supercluster comparisons (Fig. 4a,b,c, volcano plots). Several members of the Activator Protein 1 (AP1) family of transcription factors (JUN, FOS, ATF3) were consistently altered in the different superclusters. This is of importance as AP1 has been reported as one of the “core” TFs binding enhancers regulating epithelial and mesenchymal states [16, 26]. This core is subsequently modulated by other TFs such as those described in the previous paragraph. Another element strongly differing between superclusters was BHLHE40, a member of the bHLH TF family, which was found strongly downregulated in supercluster A. BHLHE40 has been reported to induce EMT as well as tumour growth and lung metastases via HBEGF exosomal release [26]. Different types of RNA-coding genes are modulated during EMT, some of which (NEAT1, MALAT1) figured in the top DEGs across the superclusters. NEAT1 was found upregulated in superclusters A and B compared to cluster C. NEAT1 has been described as promoting chemoresistance and cancer stemness [27,28,29]. More importantly, it was found to enhance glycolysis as a scaffold of key glycolytic enzymes. Its downregulation in supercluster C, which is enriched for fatty acid metabolism and the oxidative phosphorylation process, is notable. This observation may reflect a metabolic shift away from glycolysis and towards oxidative phosphorylation that was shown to exacerbate breast cancer lung micrometastases. Like NEAT1, MALAT1 is upregulated in superclusters A and B compared to C. Its effects are manifold and some controversy exists about its activities in different cancer types and settings [30, 31]. It is also interesting to note that MALAT1 has been described to interact with TEAD, which is part of the core TFs controlling the EMT process.

Finally, to relate each of the biological states to patient outcome, we selected the upregulated transcripts characteristic of each supercluster by overlapping the different contrasts (see Supplementary Table 2). The transcripts defining supercluster A correlate with better relapse-free survival (RFS) in patients suffering from basal-like breast cancer (Supplementary Fig. 4c). On the contrary, upregulated transcripts found in superclusters B and C (Supplementary Fig. 4d and e respectively) correlated with worse RFS. These results suggest that these cell states may exhibit different levels of aggressiveness driven by different sets of transcripts.