Single-cell lineage tracing of metastatic cancer reveals selection of hybrid EMT states

Abstract

The underpinnings of cancer metastasis remain poorly understood, in part due to a lack of tools for probing their emergence at high resolution. Here we present macsGESTALT, an inducible CRISPR-Cas9-based lineage recorder with highly efficient single-cell capture of both transcriptional and phylogenetic information. Applying macsGESTALT to a mouse model of metastatic pancreatic cancer, we recover ∼380,000 CRISPR target sites and reconstruct dissemination of ∼28,000 single cells across multiple metastatic sites. We find that cells occupy a continuum of epithelial-to-mesenchymal transition (EMT) states. Metastatic potential peaks in rare, late-hybrid EMT states, which are aggressively selected from a predominately epithelial ancestral pool. The gene signatures of these late-hybrid EMT states are predictive of reduced survival in both human pancreatic and lung cancer patients, highlighting their relevance to clinical disease progression. Finally, we observe evidence for in vivo propagation of S100 family gene expression across clonally distinct metastatic subpopulations.

Keywords: CRISPR; EMT; S100; barcoding; epithelial-to-mesenchymal transition; evolving barcodes; lineage tracing; metastasis; phylogenetics; single cell.

Figure 1.. macsGESTALT for high-resolution lineage tracing

(A) Genetic components of macsGESTALT. (B) Clone-level information is stored in static barcodes, while subclonal phylogenetic information is dynamically encoded into evolving barcodes via insertions and deletions (indels, blue and red bars) induced by doxycycline. (C) Two example clones from a population with n clones, each with a random number of integrated barcodes. Evolving barcode edits are encoded and inherited as cells divide. (D) Generation of a macsGESTALT barcoded population of cells and experimental workflow. (E) macsGESTALT analysis workflow. Dox, doxycycline; rtTA, reverse tetracycline transactivator; TRE, tetracycline-responsive element. See also Figures S1 and S2.

Figure 2.. Most metastases arise from rare, transcriptionally distinct clones

(A) Schematic of metastasis lineage tracing model. (B) Clonal reconstruction using static barcodes, where clones are numbered by size in the primary tumor. Percentage contribution to each harvest site (circle size) and enrichment compared with the primary tumor (circle color) are visualized. Top annotations show each clone’s Leiden transcriptional cluster and aggression assignments as in (H) and (I), respectively. (C) Cumulative fraction of each clone in each disseminated site (red) and primary tumor (black). Dotted lines represent the theoretical scenario of perfect clone size equality. (D) UMAP plot of 28,028 single cells containing both lineage and transcriptional information. Cells are colored by clone, with select large clones highlighted (as mouse.clone). (E and F) Two representative non-aggressive clones. (E) M1.13, (F) M2.10. (G) A representative clone of medium aggression. (H) Leiden transcriptional clustering of (D). (I) Cells colored by clonal aggression. (J) Number of non-, mid-, or high-aggression clones of 95 total. See also Figures S3 and S4.

Figure 3.. A transcriptional EMT continuum in vivo

(A) UMAP plot of M1, colored by clone, with the five largest clones annotated. Circled region indicates the transcriptional space where smaller, non-aggressive clones reside. (B–G) Expression of canonical epithelial (B, Epcam; C, Muc1; D, Cdh1) and mesenchymal (E, Sparc; F, Zeb2; G, Col3a1) markers. (H) Unbiased trajectory inference revealing a pseudotime axis matching EMT (pseudoEMT). (I) Expression of (B–G) plotted along pseudoEMT and colored by clone as in (A). (J) Hierarchical clustering of kinetic curves for the top 3,000 differentially expressed genes across pseudoEMT (q = 0, Moran’s I > 0.1). Gene clusters are labeled from epithelial, E, to hybrid, H1–H4, to mesenchymal, M, based on expression across pseudoEMT. Gene set analysis using MSigDB Hallmarks for each gene cluster (hypergeometric test, p < 0.05). Oxphos, oxidative phosphorylation. (K) Significantly enriched motifs (hypergeometric test, p < 0.05) in promoters for each gene cluster, with canonical EMT master regulators highlighted. See also Figure S5 and Tables S1–S3.

Figure 4.. High-resolution subclonal lineage reconstruction of metastatic cancer

(A) Percentage at which each base is mutated in 76,974 evolving barcodes across both mice. Target-site spacers (light gray) and PAMs (dark gray). (B) Edit types observed at each target site. (C) Example phylogenetic reconstruction of a small clade within clone M1.1. Clade M1.1.310 (root node in red) contains six distinct subclones composed of 58 cells from five different harvest sites. Each cell in this clade has six evolving barcodes, illustrated by white bars with edits colored as in (B). Cells with the same barcode editing pattern are grouped into a subclone (terminal black nodes) and dissemination (E_H) is quantified. For each subclone, individual cells are stacked and colored by their harvest site on the far right. (D) Circle packing plot of the full single-cell phylogeny of M1, with clade M1.1.310 from (C) circled in red. Outermost circles define clones, with the first six clones labeled. Within each clone, nested circles group increasingly related cells. Innermost circles contain cells from reconstructed subclones. Each point represents a single cell, colored by harvest site. (E) Cumulative fraction of each subclone of clone M1.1 in each harvest site. Dotted line represents perfect subclone-size equality. See also Figure S3.

Figure 5.. Peak metastatic aggression corresponds to late-hybrid EMT states

(A and B) Circle packing plots of the phylogenetic structure of clone M1.1 with subclones colored by mean pseudoEMT (A) and by dissemination score (B). (C) Relationship between metastatic dissemination and pseudoEMT for subclones from (A and B). (D) Density along pseudoEMT of M1.1 cells and their increasingly ancestral (arrow) phylogenetic groupings, examples of which are highlighted in (A). (E) Relationship between PDAC patient survival (TCGA-PAAD, n = 173) and patient enrichment scores for each pseudoEMT gene cluster using Cox regression analysis, with the hazard ratio for each gene cluster displayed (*p < 0.05, •p < 0.1). Square sizes are inversely proportional to p value. See also Table S4.

Figure 6.. A process complementary to canonical EMT

(A) Lineage tree for M2 subclones, where branches and nodes are colored by clone and scaled by the number of cells they relate. (B) UMAP of M2 cells, colored as in (A), with five large, aggressive clones labeled, as well as M2.1 (green), which was the largest clone in the primary tumor but poorly metastatic. Circled region indicates the transcriptional space where smaller, non-aggressive clones reside. (C) Relationship between PDAC patient survival (TCGA-PAAD, n = 173) and enrichment scores for genes associated with subclonal dissemination using Cox regression analysis (**p < 0.01), with the hazard ratio displayed. Square sizes are inversely proportional to p value. (D–H) Canonical epithelial (D, Ocln; E, Epcam; F, Lgals4) and mesenchymal (G, Prrx1; H, Zeb2) markers. (I and J) Markers with inconsistent expression patterns in the dominant clone, M2.2 (I, Sparc; J, Muc1). (K) Highly expressed genes ranked by association (q < 0.05) with subclonal dissemination. (L) Aggregated single-cell gene expression of the S100a family for each clone, colored by aggression (as defined in Figure 2B) and grouped by mouse. Intramouse comparisons between dominant/aggressive clones versus all others are indicated above each violin. Comparisons between mice for all clones (black) and only dominant/aggressive clones (red) are indicated above the line (Welch’s t test, ****p < 0.0001, ***p < 0.001, ns, not significant). See also Figure S6 and Table S5.