Progressive plasticity during colorectal cancer metastasis

Abstract

As cancers progress, they become increasingly aggressive-metastatic tumours are less responsive to first-line therapies than primary tumours, they acquire resistance to successive therapies and eventually cause death^1,2. Mutations are largely conserved between primary and metastatic tumours from the same patients, suggesting that non-genetic phenotypic plasticity has a major role in cancer progression and therapy resistance^3-5. However, we lack an understanding of metastatic cell states and the mechanisms by which they transition. Here, in a cohort of biospecimen trios from same-patient normal colon, primary and metastatic colorectal cancer, we show that, although primary tumours largely adopt LGR5⁺ intestinal stem-like states, metastases display progressive plasticity. Cancer cells lose intestinal cell identities and reprogram into a highly conserved fetal progenitor state before undergoing non-canonical differentiation into divergent squamous and neuroendocrine-like states, a process that is exacerbated in metastasis and by chemotherapy and is associated with poor patient survival. Using matched patient-derived organoids, we demonstrate that metastatic cells exhibit greater cell-autonomous multilineage differentiation potential in response to microenvironment cues compared with their intestinal lineage-restricted primary tumour counterparts. We identify PROX1 as a repressor of non-intestinal lineage in the fetal progenitor state, and show that downregulation of PROX1 licenses non-canonical reprogramming.

Fig. 1. Non-canonical transcriptional programs in CRC are associated with metastasis and poor outcomes.

a, Study design. Matched biospecimen trios of normal colon, primary CRC and metastasis were collected from 31 patients, processed fresh for single-cell transcriptomics and organoid generation, and formalin-fixed and paraffin-embedded for multiplexed immunofluorescence analysis. Organoids were used for functional studies in vitro or through orthotopic xenotransplantation into the caecum or liver. b, Hotspot modules in all CRC tumour cells. The heat map comprises 2,003 highly variable genes with significant autocorrelation (FDR < 0.01), grouped into 4 canonical and 6 non-canonical CRC-derived gene modules (Methods and Supplementary Table 4). Dev., development. c, The distribution of module proportions in metastatic tumours from treated (red dots) and untreated (green dots) patients; module labels are based on the >0.75 quantile score in a given cell for the gene module (Methods). The vertical line divides canonical intestinal (right) from other (left) cell types. d,e, The log-ratio of metastasis- to primary-derived tumour module proportions in each patient sample, based on the accumulation of non-canonical (d) or canonical (e) modules. Metastatic tumours are significantly enriched for cells expressing non-canonical modules, whereas primary tumours are enriched for cells expressing canonical modules (P = 0.001, one-sided rank-sum test; Methods). f,g, Associations of module enrichments with tumour recurrence (f) and the survival status (g) of patient donors in two independent clinical cohorts. Associations are shown for 108 patients with rectal adenocarcinoma (LARC) and for 445 patients with colon adenocarcinoma (TCGA). Enrichment scores were calculated using ssGSEA of bulk transcriptomic data from each patient (Methods). Each gene module is plotted according to the Mann–Whitney U-test statistic of distal recurrence (x axis) and new tumour event or survival under 10 years (y axis). Abs., absorptive intestine; endo., endoderm development; inj., injury repair; int., intestine; neuro., neuroendocrine; osteo., osteoblast; sec., secretory intestine; squa., squamous.

Fig. 2. A conserved fetal progenitor intermediate bridges canonical and non-canonical states in CRC.

a, Force-directed layout of cells from liver metastasis (met.) from patient KG146 (1,279 cells), showing a diversity of canonical and non-canonical states. Each cell is coloured according to its maximum module score. b, Trends in gene module scores along DC1 observed in all tumour cells from patient KG146 (top) (Methods). Each row depicts the module score along DC1 from the 20th percentile value (white) to the maximum value (highest saturation). Expression of the fetal signature peaks before the non-canonical peak. The closed and open arrowheads correspond to the 75th percentile of fetal and predominant terminal non-canonical module scores, respectively. Bottom, positions of tumour cells along DC1. c, Uniform manifold approximation and projection (UMAP) embedding of human fetal colon cells. First-trimester (6–11 weeks after conception) and second-trimester (12–17 weeks after conception) cells are coloured by cell type according to published annotations. TA, transit-amplifying cells. d–f, The faction of cells expressing >0.75 quantile score for a given module (or the core fetal signature, maroon rectangle) in samples from the four patients with most non-canonical cells (KG146, KG182, KG150 and KG183); d), all samples in our cohort (e) and samples from the ref. cohort (f). Shared module expression reveals a consistent progression from canonical to non-canonical fates across patients and cohorts.

Fig. 3. Organoid models reveal distinct contributions of tumour and microenvironment to metastatic plasticity.

a, Bright-field microscopy showing the morphology of paired primary rectal-tumour-derived (OKG146P) and liver-metastasis-derived (OKG146Li) organoids grown in HISC medium, 7 days after seeding as single cells (2,000 cells per 40 µl Matrigel). Scale bars, 500 μm. b, Cell-state assignment probabilities per cell, calculated as Markov absorption probabilities (Methods) for OKG146P and OKG146Li organoids grown in HISC (left) and IGFF (right) medium. The lines indicate density contours within the 5th and 95th percentiles of each distribution, and the dots indicate individual cells. Diff. int., differentiated intestinal states. c,d, The normalized average radiance was measured by weekly ex vivo bioluminescence imaging after caecal injection of 200,000 cells (c) or intrahepatic injection of 500,000 cells (d) from OKG146P (primary tumour) and OKG146Li (metastasis) lines in NSG mice, normalized to the signal immediately after injection (week 0). Data are mean ± s.e.m. of n = 6 (OKG146P) and n = 7 (OKG146Li) mice (c) and n = 5 (OKG146P) and n = 4 (OKG146Li) mice (d). For c and d, statistical analysis was performed using two-sided Mann–Whitney rank-sum tests; P = 0.2246 (c) and 0.0143 (d) comparing between signal at end point. Source Data

Fig. 4. PROX1 encodes a fetal-state-associated transcription factor that inhibits non-canonical transdifferentiation.

a, The probability of classifying cells as non-canonical in scRNA-seq data from 146 primary (left) and liver metastasis (right) shControl and shPROX1 organoid lines, 7 days after induction with doxycycline. Cells were classified as canonical or non-canonical using a manifold-based classifier that combines methods from Harmony and PhenoGraph (Methods). The vertical black line indicates the median. Statistical analysis was performed using two-sided t-tests; ***P < 0.001; P = 0.0 (left) and P = 2.31 × 10⁻⁷ (right). b, The relative expression of canonical and non-canonical differentiation markers in organoid lines expressing shRNAs targeting PROX1 or control shRNA (shCtrl) cultured in HISC medium containing 2 μg ml⁻¹ doxycycline for 7 days. Quantitative PCR with reverse transcription (RT–qPCR) data are normalized to GAPDH mRNA expression. n = 4 replicates per group. Statistical analysis was performed using two-sided t-tests with Benjamini–Hochberg correction; *P < 0.05. c, Model of cell-state transitions during metastatic tumour progression in CRC. Cancer cells in the primary tumour first enter an ISC-like state, then cells at the tumour invasion front undergo developmental reversion into a fetal progenitor-like state, enabling differentiation into divergent non-canonical states, including neuroendocrine- and squamous-like, that are enriched during metastatic outgrowth. Entry into the highly plastic fetal progenitor state is triggered by epithelial injury during tumour dissemination or after therapy, allowing tumour regenerative cells to express non-canonical gene programs and adapt to diverse stresses. Induction of PROX1 inhibits non-canonical gene expression in injured normal epithelia, enabling tissue regeneration, whereas PROX1-responsive intestinal lineage restriction is progressively lost during cancer progression, licensing non-canonical differentiation. Source Data

Extended Data Fig. 1. Clinical and genomic characteristics of metastatic colorectal cancer biospecimen cohort.

a, Clinical characteristics of patients in the study cohort, including key demographic, clinical, metastatic site, treatment and outcome variables, as well as whether scRNA-seq data was collected from each biospecimen or organoid (also see Supplementary Table 1). Longitudinal samples refer to additional metachronous tumour samples collected at time of progression, subsequent to initial synchronous tissue collection. No more than one biospecimen was sequenced per site. b, Genomic features of tumours sequenced using the MSK-IMPACT platform. c, Summary of biospecimens and datasets collected.

Extended Data Fig. 2. Clinical and genomic features of patient cancers and collected biospecimens.

a, Expression of canonical marker genes (rows) across all cells (columns) grouped by compartment, cell type and cluster. Colour values represent normalized transcript counts (Methods). b, Uniform manifold approximation and projection (UMAP) embedding of all cells coloured by immune (111,609 cells), epithelial (47,437 cells) or stromal (5,258 cells) compartment. c, Copy number changes calculated with InferCNV and binned by genomic region, per epithelial cell cluster (cells clustered based on inferred copy number matrix; Methods). Rows represent mean inferred copy number relative to diploid reference population of non-tumour epithelium, for each cluster of cells. Cancer cells are called on a per-cluster basis according to their mean copy number profile; right column, green (no CNV) or red (CNV). d, Inferred copy number of tumour cells for a representative patient (top, mapping to three inferCNV clusters, including a clearly dominant cluster with 82.3% of cells) compared to copy number values estimated by FACETS using bulk sequencing for a targeted gene panel (bottom; see Methods). e, UMAP of all epithelial cells as in Fig. 1d, coloured by CNV classification. f, UMAP of all non-tumour epithelial cells (21,297 cells), coloured by cell type annotation (Methods). g, Expression of colon epithelial cell type markers across annotated cell types. Dot size scales with proportion of cells in a cell type that express each gene; colour intensity indicates mean z-scored, log-normalized expression of each gene. Genes (rows) are coloured by the cell type (as in f) they mark. h, UMAP of all epithelial cell scRNA-seq profiles collected from non-tumour colon (20,817 cells), CRC primary tumour (14,402 cells), and CRC metastatic tumour (12,218 cells) samples, coloured by cell type. ISC, Intestinal stem cell. i, UMAP of epithelial cells, coloured by tissue origin. j, UMAP of epithelial cells, coloured by pre-surgical treatment status (see Supplementary Fig. 1 for patient treatment details).

Extended Data Fig. 3. Untreated CRC tumours express early stem cell and mixed lineage programs.

a, Force-directed layout (FDL) of all treatment-naive epithelial cells (13,935 cells), coloured by tumour type (primary, metastasis) for tumour cells or cell type for all other cells. b, Major cell types of the colon crypt. c, FDL of treatment-naive epithelial cells as in a, coloured by our de novo ISC gene signature score (Methods). Scale indicates average, z-normalized gene expression (Methods). d, Distributions of first principal component (PC1) values (top) or relative ISC signature expression (bottom) for all treatment-naive epithelial cell types (13,935 cells) ordered by median PC1 value and coloured as in a and b. Gene set enrichment analysis (GSEA) performed on genes ranked by their PC1 loadings was used to determine ‘differentiated’ versus ‘stem-like’ directions (Methods; see Supplementary Table 2 for all gene sets and enrichment scores). Dots, median values; lines, interquartile range. e, Mean expression of ISC gene signature genes, z-normalized per gene across cells, within all treatment-naive epithelial cell types (Supplementary Table 2 and Methods). Rows organized by cell type and columns organized by clustering on the genes. f, Genes differentially expressed between ISCs and all treatment-naive tumour cells. Top-enriched gene sets according to GSEA performed on genes ranked by differential expression in tumour cells (Methods); highlighted genes are from leading edge subsets. OXPHOS, oxidative phosphorylation; TA, transit amplifying; ISC, intestinal stem cell; WNT, WNT-beta-catenin signalling; EMT, epithelial-mesenchymal transition; UPR, unfolded protein response; MTORC1, MTORC1 signalling. logFC, log fold-change. -log10 Adj. P-Val., two-sided t-test, –log₁₀-adjusted p value with Benjamini-Hochberg correction. g, Gene-gene correlations in cells of normal intestine (upper diagonal), and untreated tumours (lower diagonal) showing expression of both ISC and mixed lineage markers in treatment-naive tumour cells. Genes (some highlighted) consist of the top 100 differentially expressed genes (DEGs) in ISCs, enterocytes, and goblet cells compared to all other non-tumour epithelial cells. Colour values represent average gene-gene correlations calculated for each sample (Methods). h, Abundance of ISC- and enterocyte-specific DEGs in normal ISCs and enterocytes, and primary tumour and metastasis cells from untreated patients (Methods). Density plot only includes markers with top quartile expression per gene across all tumour cells. i, Same as h, for enterocyte and goblet cell markers. j, Same as i, for ISC and goblet cell markers.

Extended Data Fig. 4. Decreased expression of ISC programs in CRC metastasis.

a, Histograms of ISC gene signature scores for all tumour cells (26,145 cells). Dots, median; interquartile range (IQR), black bars; 1.5x IQR, black lines. ***, untreated p = 9.75e-15; ***, treated p = 1.11e-157 two-sided Wilcoxon rank-sum test. b, Vectra multiplex immunofluorescence from matched trios of normal colon, primary tumour and liver metastasis, showing OLFM4 (ISC marker), Ki67 (proliferation marker) and DAPI in representative samples from previously untreated (green) and treated (brown) patients. c, Histograms quantifying OLFM4-high expressing cells (Methods; untreated primary = 3,61,532, treated primary = 1,180,237, untreated metastasis = 565,014, treated metastasis = 150,1084). Dots, median; IQR, black bars; black lines, 1.5x IQR; *, primary p = 0.042; ***, treated p = 0.003; ***, metastasis p = 9.50e-9; two-sided Wilcoxon rank-sum test.

Extended Data Fig. 5. Hotspot identifies canonical intestinal and non-canonical gene expression modules in tumour cells from the full metastatic CRC cohort.

a, Average pairwise local correlation score (Methods) between Hotspot modules before grouping. Each entry is the average local correlation between all pairs of genes assigned to the corresponding module pair (row and column). Numbered squares are coloured according to final module annotation after grouping. b, z-normalized gene expression within the top-scoring cells for each ungrouped Hotspot module (rows) (Methods). Columns are coloured according to the module assigned to the gene on the x-axis. Dot size indicates proportion of cells in a cluster that score for that gene module; colour intensity indicates mean module score of that cluster. Numbered squares are coloured as in a. c, z-normalized gene expression within groups of top-scoring cells for each grouped Hotspot module as in b. d, Kernel density plots depicting entropy across patients within high-scoring tumour cells for each module (Methods). High entropy indicates that cells with high score for a module (>1 s.d. above mean, calculated across all tumour cells) come from a diverse set of patient samples; low entropy indicates that cells primarily originate from a single patient. e–j, FDL of all tumour cells (26,145 cells). Each cell is coloured by gene module score for 4 canonical (e) or 6 non-canonical (f) CRC gene modules, tumour site (g), treatment status (h), according to its maximum module score, or grey if no score exceeds the module’s 50th percentile across all cells, coloured by the number of Hotspot modules for which it scores higher than the 0.5 (i) or 0.75 (j) quantile module score across all cells.

Extended Data Fig. 6. Shift from canonical intestinal marker expression in metastases and therapy-treated samples is associated with poor clinical outcomes.

a,b, Vectra multiplex immunofluorescence from additional representative matched trios stained for canonical intestinal and non-canonical markers, a) OLFM4 (ISC marker), CK5 (squamous cell marker), TROP2 (injury repair state marker), Ki67 (proliferation marker) and DAPI (nuclei), and b) HER2 (pan-epithelial), CDX2 (intestinal lineage), CK20 (differentiated intestinal), SOX2 (multipotential stem cell) and CHGA (neuroendocrine marker), demonstrating increasingly disordered morphology and increasing expression of non-canonical markers during tumour progression. Patient identifiers (columns) indicate treated (TR, red) or untreated (UT, green) status. c–f, Fraction of tumour cells per field of view with high expression (‘high expression’ for each marker based on minimal expression threshold determined by knee-point decile value from all marker-positive cells) of c) CK20 (differentiation; *, treatment-naïve p = 0.025; ***, treated p = 2.92e-4; *, primary p = 0.0286; **, metastasis p = 0.002), d) OLFM4 (intestinal stem cell; ***, treated p = 3.24e-4; *, primary p = 0.042; ***, metastasis p = 8.50e-9), e) TROP2 (injury repair; ***, treatment-naïve p = 4.65e-15; *, treated p = 0.017; ***, primary p = 1.06e-9) and f) CK5 (squamous; **, primary p = 0.002; ***, metastasis p = 7.80e-6). Orange, primary tumours; purple, metastases; white line, median; black bars, IQR; black lines, 1.5x IQR. n = 2,096,785 cells (c), n = 1,408,818 cells (d–f), rank-sum test. g, Cumulative fraction of cells expressing >0.75 quantile score for a given module in 5 patients from ref. (p = 0.008). One sided rank-sum test (Methods). h,i, Enrichment of gene module scores in scRNA-seq data with respect to h) patient baseline clinical, pathological and treatment attributes and i), CRC consensus molecular subtype (CMS) classifications; Mann–Whitney U rank-sum test. j, ssGSEA gene module score enrichment with respect to patient baseline clinical, pathological and treatment attributes for 445 patients with colon (COAD) adenocarcinoma from the TCGA cohort; Mann–Whitney U test. k, Kaplan-Meier plots showing disease-free survival for patients in the TCGA cohort with high or low ssGSEA enrichment for the indicated signatures. A patient is signature-high if the enrichment score of their tumour is >1 s.d. above the mean, calculated across all patients, and signature-low if it is <1 s.d. below the mean (log-rank; n = 147, 158, 144, and 123 for Absorptive p = 0.003, Secretory p = 0.003, Endoderm p = 0.001, and Squamous p = 0.026, respectively).

Extended Data Fig. 7. Cell state progression is conserved across patients.

a, Top, trends in gene module scores along DC1, observed in all tumour cells from patients KG152, KG180, KG183, and patient s1321 from ref. (Methods). Rows depict module score along DC1 from 20th percentile (white) to maximum value (highest saturation). Fetal signature expression peaks between canonical and NC trends; the closed and open arrowheads correspond to the 75th percentile of fetal and predominant terminal non-canonical module scores, respectively. Bottom, positions of tumour cells along DC1. b, Jaccard index (similarity metric defined as no. shared genes / total no. genes in both sets) between the fetal signature and published signatures from mouse (left), and normalized overlap with the 14 shared genes from (d) (right). c, Average z-normalized gene set scores in the cell types of the human fetal dataset, demonstrating overlap of our data with the human fetal signature, compared to previous mouse signatures. d, Overlap of fetal signature genes with Pearson correlation >0.5 across the four patient samples harbouring a substantial number of cells in non-canonical differentiated states. e, Relative correlation of core fetal signature with different modules, showing highest correlation for Endoderm Development module across all tumours. f, Log-ratio of metastasis- to primary-derived tumour cells that exhibit >0.75 quantile score for the 14-gene conserved fetal signature, for each patient sample (p = 0.0004, one-sided rank sum test; Methods). g, Kaplan-Meier plots showing disease-free survival for patients in the LARC cohort with high or low ssGSEA enrichment for the fetal progenitor signature. A patient is fetal-progenitor-high if their enrichment score is >1 s.d. above the mean, calculated across all patients, and fetal-progenitor-low if it is <1 s.d. below the mean (log-rank, n = 42, p = 0.0350). h, Same as g, for patients in the TCGA cohort (log-rank, n = 134, p = 0.0016).

Extended Data Fig. 8. Gene expression changes along Palantir pseudotime for squamous and neuroendocrine trajectories.

a, UMAP of 3351 cells from KG146 primary tumour, synchronous liver metastasis, and metachronous lung metastasis. Cell are coloured by maximum module score, or grey if no module score exceeds its 25th percentile (top-left); by Palantir pseudotime (top-middle); by Palantir branch probability for the squamous-annotated (top-right) and neuroendocrine-annotated (bottom-left) branches; or by the fetal signature (bottom-left). b, Pearson correlations of Palantir branch probabilities with the expression of genes associated with squamous (top) and neuroendocrine (bottom) transdifferentiation, across tumours from patients KG146, KG182 and KG150. c, Gene expression trends along Palantir pseudotime for selected signalling pathways during trans-differentiation of KG146 tumour cells toward squamous or neuroendocrine-like terminal states. Trends are computed using generalized additive models (GAMs) in Palantir. Solid lines represent mean expression, shaded regions represent 1 s.d.

Extended Data Fig. 9. Transcriptomic comparison of patient tumours and organoids deconvolves cell-autonomous and microenvironmental contributions to cell state.

a, MSK-IMPACT results showing concordant mutation profiles of primary tumour and patient-derived organoids for patient KG146. b, FDL of KG146 primary CRC tumour cells (n = 880) coloured by cluster-level cell state annotations (Methods). c, Same as b, for KG146Li patient tumour cells (n = 1,279). d, CRC gene module scores within cell clusters for patient tumour and organoid scRNA-seq samples (Methods). Dot size indicates proportion of cells in a cluster that score for that gene module; colour intensity indicates mean module score of that cluster. Cluster assignments are made within each sample, and scores are calculated separately for all patient tumour samples and all organoid samples. Squares above patient tumour samples indicate the cluster classification (blue, differentiated intestine-like states; green, ISC-like and TA-like states; red, non-canonical states). e,f, Relative expression of differentiated intestine, ISC, fetal and NC differentiated states in patient-derived organoid line OKG182CW2 (secondary metastasis from chest wall) (e) or OKG183Li2 (secondary metastasis from liver) (f) cultured for 7 d in IGFF compared to HISC media. RT-qPCR normalized to GAPDH mRNA expression level (n = 4 technical replicates, Mann–Whitney rank-sum test, *, p < 0.05, **, p < 0.01, ***, p < 0.001****, p < 0.0001; KG182CW2 p-values, KRT20: 8.09e-09, FABP1: 6.82e-07, TFF3: 0.001, LGR5: 2.51e-14, GPC1: 0.003, TMEM132A: 1.36e-06, NEUROD1: 0.02, CHGB: 3.05e-06, KRT23: 1.24e-06; KG183Li2 p-values, KRT20: 1.59e-4, FABP1: 0.012, TFF3: 5.120e-10, LGR5: 1.22e-10, GPC1: 1.82e-4, TMEM132A: 5.43e-06, KRT23: 1.22e-09. g,h, Normalized average radiance measured by weekly ex vivo bioluminescence imaging following intrahepatic injection of 500,000 cells in NSG mice, normalized to signal immediately following injection (week 0), for g) OKG146P (primary tumour) and OKG146Li (metastasis) lines (n = 5 and 4, respectively; p = 0.0143, Mann–Whitney two-sided rank-sum test at week 12; p = 1.458e-13, normalized radiance > 10, Fisher’s exact test, all weeks), and h) OKG136P (primary tumour) and OKG136Li (metastasis) lines (n = 10 and 10; p = 0.39, Mann–Whitney two-sided rank-sum test at week 12; p = 2.35e-6, normalized radiance > 10, Fisher’s exact test, all weeks). Error bars, s.e.m. i Representative Vectra images of OKG146Li intrahepatic xenografts, harvested at 13 weeks post injection and stained for (left) DAPI (nuclei), VIM (vimentin, stroma), HER2 (epithelial differentiation), CDX2 (intestine), CK20 (differentiated intestine), SOX2 (stem cell) and CHGA (neuroendocrine) markers, and for (right) DAPI (nuclei), OLFM4 (ISC), CK5 (squamous), Ki67 (proliferation) and TROP2 (injury repair state) markers. Channels are shown for all markers (left), intestinal markers (centre) and SOX2 and CHGA (right); 20x magnification corresponds to 10x inset. Source Data

Extended Data Fig. 10. Chemotherapy induces expression of non-canonical programs.

a, Gene module scores within cell clusters for patient tumour and organoids cultured in HISC media and treated with irinotecan, as in Extended Data Fig. 9d. b, Cell state assignment probabilities per cell for live cells from OKG146 organoids grown in HISC media and treated with 250 nM irinotecan chemotherapy for 7 d. Lines indicate density contours within 5th and 95th percentiles of each distribution, dots indicate individual cells. c, Fetal signature score distributions in primary (orange) and metastatic (purple) cancer cells from patient KG146 in patient tumour (left) and patient-derived organoids cultured in IGFF, HISC or HISC media containing 250 nM irinotecan. Fetal signature is the same as in Fig. 4 and calculated on z-normalized, imputed expression. Dots, median; IQR, black bars; black lines, 1.5x IQR; ***, p < 0.001, two-sided t-tests comparing primary and metastasis samples between conditions. p = 3.48e-26, 2.82e-06, 7.40e-147 for patient, IGFF: HISC respectively (n = 12,016 cells). d, Relative expression (log fold-change) of NC differentiation markers in primary tumour (OKG136P, OKG173P, OKG146P) and metastasis (OKG136Li, OKG173Li, OKG146Li) organoids cultured for 7 d in HISC media supplemented with 250 nM irinotecan compared to HISC media alone. RT-qPCR normalized to GAPDH expression (n = 4 technical replicates, two-sided t-test; *, p < 0.05). Source Data

Extended Data Fig. 11. Identification of PROX1 as a fetal-state associated transcription factor.

a, Gene trends along aligned canonical–non-canonical axes, of the 6 top-ranked fetal-associated human TFs in four patients. TFs are ranked according to their mean correlations of imputed gene expression with the fetal signature score in these patients. b, Relative expression, from scRNA-seq data, of top 6 fetal-associated TFs in OKG146P (light bars) and OKG146Li (dark bars) organoids cultured for 7 d in HISC media containing 250 nM irinotecan, compared to HISC media alone. c, Pearson correlation of PROX1 expression with module gene expression across this cohort or the ref. cohort. Genes are grouped by module and ordered within each module by hierarchical clustering; modules are ordered by average correlation with PROX1 expression. d, Hematoxylin and eosin (H&E) and immunofluorescence imaging of PROX1, CDX2 (intestinal lineage) and VIM (vimentin, stroma), showing PROX1 expression in poorly differentiated cells along the invasion fronts of two patient primary tumors.

Extended Data Fig. 12. PROX1 is a context-dependent stabilizer of intestinal identity in an epithelial injury-associated, fetal-like state.

a, Representative images and b) quantification of whole mount immunofluorescence of OKG146Li organoids expressing shRNAs targeting PROX1 or control, cultured in HISC media with 2 μg/ml doxycycline for 7 d, stained for PROX1 (red) and DAPI (blue). Scale bar, 100 μm. Boxes, IQR; whiskers, 1.5 * IQR. n = 10, 12 fields of view for shCtrl; 12, 12 for shPROX1-4; 12, 13 for shPROX1-5, quantified per HISC or HISC + 2 μg/ml doxycycline condition, respectively. ****, p < 0.0001, ns = not significant. shControl, p = 0.79, shPROX1-1, p = 5.62e-07, shPROX1-2, p = 1.20e-06, two-sided t-test. c, PROX1 protein expression in OK146P organoids expressing doxycycline-inducible shRNAs targeting PROX1 or control, cultured in IGFF medium with 2 μg/ml doxycycline for 7 d. GAPDH, loading control. See Supplementary Fig. 4 for original gel image. d,e, Normalized average radiance measured by weekly ex vivo bioluminescence imaging following intrahepatic injection of 500,000 organoids from OKG146P (d) or OKG136P (e) lines bearing shCtrl, shPROX1-1 and shPROX1-2 (n = 9, 10 and 10, respectively, for each organoid line from each patient) in NSG mice, normalized to signal immediately following injection (week 0). Error bars, s.e.m. ** p = 0.008, ***, p = 0.001, two-sided Mann–Whitney rank-sum test. p = 0.0069 for OKG146P shPROX1-1 and p > 0.05 for all others, normalized radiance > 10, Fisher’s exact test, all weeks. f, Representative images of OKG146P and OKG146Li organoids expressing doxycycline-inducible shRNAs targeting PROX1 or control shRNA, 7 d after plating 750 single cells/15 μl Matrigel in HISC + 2 μg/ml doxycycline (scale bar = 1000 μm). g, Organoid initiation capacity (number of organoids formed per 750 single cells/15 μl Matrigel) of organoids in f). n = 8 independent replicates. Boxes, IQR; whiskers, 1.5 * IQR. *, p < 0.05, ****, p < 0.0001. shPROX1-1 p = 0.26, 2.9e-08, 0.010, 0.00018, 0.68, 0.93, 0.90; shPROX1-2 p = 9.1e-05, 0.018, 0.0023, 0.0025, 0.069, 0.073, 0.02 for: KG136P, KG136Li, KG173Li, KG146P, KG146Li, KG182CW2, KG183Li2. One-sided t-test, Benjamini-Hochberg correction. h, Relative expression, from scRNA-seq data, of the conserved fetal signature in OKG146P (n = 3098), OKG146Li shControl (n = 2840) and OKG146P (n = 3504), OKG146Li (n = 2794) shPROX1 organoid lines, 7 d after induction with dox. Dots, median; IQR, black bars; black lines, 1.5x IQR, p = 0.054, two-sided t-test. Source Data