Oncogenic enhancers prime quiescent metastatic cells to escape NK immune surveillance by eliciting transcriptional memory

Abstract

Metastasis arises from disseminated tumour cells (DTCs) that are characterized by intrinsic phenotypic plasticity and the capability of seeding to secondary organs. DTCs can remain latent for years before giving rise to symptomatic overt metastasis. In this context, DTCs fluctuate between a quiescent and proliferative state in response to systemic and microenvironmental signals including immune-mediated surveillance. Despite its relevance, how intrinsic mechanisms sustain DTCs plasticity has not been addressed. By interrogating the epigenetic state of metastatic cells, we find that tumour progression is coupled with the activation of oncogenic enhancers that are organized in variable interconnected chromatin domains. This spatial chromatin context leads to the activation of a robust transcriptional response upon repeated exposure to retinoic acid (RA). We show that this adaptive mechanism sustains the quiescence of DTCs through the activation of the master regulator SOX9. Finally, we determine that RA-stimulated transcriptional memory increases the fitness of metastatic cells by supporting the escape of quiescent DTCs from NK-mediated immune surveillance. Overall, these findings highlight the contribution of oncogenic enhancers in establishing transcriptional memories as an adaptive mechanism to reinforce cancer dormancy and immune escape, thus amenable for therapeutic intervention.

Fig. 1. Recapitulating metastatic progression of TNBC ex-vivo.

a Schematic representation of model derivation and analysis approach (created with BioRender.com). b Representative images of IHC and quantification of metastatic size (DTCs: n < 10 cells; micro-metastasis: 10 <  n < 100 cells; macro-metastasis: n > 100 cells), p27 and Ki67 positive cells in lung and liver metastases, and primary tumor. Quantification was performed on n = 4–16 Field of Views (FOVs), merged from three independent experiments. For liver and lung tissues, a representative image showing both micro- and macro-metastases is shown, with corresponding zoomed-in crops (black: micro-metastasis; gray: macro-metastasis). Red scale bar = 1 mm; blue scale bar = 100 µm; black scale bar = 40 µm. c Representative images of tIMEC, XD, and MD tumoroids embedded in collagen and quantification of seeding cells and distance from the spheroid. Calcein AM (green); H2B-mCherry (red); scale bar: 100 µm. d Representative images of tIMEC, XD, and MD tumoroids embedded in collagen and quantification of Col-¾ mean intensity. Calcein AM (green); H2B-mCherry (red); Col-¾ (purple) scale bar: 100 µm. e Representative images and quantification of Ki67 positive cells in tIMEC, XD, and MD cells. Ki67, white; DAPI, blue; scale bar = 10 µm. f Representative images and quantification of p27 positive cells in tIMEC, XD, and MD cells. p27, white; DAPI, blue; scale bar = 10 µm. g Representative histograms of FACS analysis on dye-retaining (APC-A) cells and quantification of APC-A+ cells. h Representative images and quantification of SOX9 nuclear mean intensity in tIMEC, XD, and MD cells. SOX9, white; DAPI, blue; scale bar = 10 µm. i Representative images of IHC and quantification of SOX9 positive cells in primary tumor, lung, and liver metastases and in subgroups (DTCs, micro- and macro-metastasis) of lung metastasis. Quantification was performed on n = 4–24 Field of Views (FOVs), merged from three independent experiments. The barplots in (b), (f–i) are means of three independent biological replicates ± S.E.M. The barplots in (e) are means of four independent biological replicates ± S.E.M. The box plots in (d) indicate median values (middle lines), first and third quartiles (box edges), and 10th and 90th percentiles (error bars) retrieved from three independent biological replicates. The violin plots in (c) and (h) indicate median values (middle lines), first and third quartiles (dashed lines) retrieved from three independent biological replicates. Statistical significance was determined by one-tailed unpaired student’s t-test.

Fig. 2. Metastatic onset is associated with epigenomic rewiring.

a Multidimensional scaling (MDS) plots of batch-corrected logCPM values with samples (tIMEC, XD, PT, and MD) colored by cell type. Distances on the plot correspond to the leading fold change, which is the average (root-mean-square) log2fold change for the 5000 most divergent peaks between each pair of samples. b Uniform Manifold Approximation and Projection (UMAP) projections of tIMEC, XD, and MD accessibility data along with publicly available accessibility datasets representing human primary cancer samples from different locations, human primary non-cancerous tissues from different locations. A zoomed-in crop highlighting tIMEC, XD, and MD samples is shown. c Boxplot showing Spearman pairwise distance between tIMEC, XD, MD, IMEC samples, and the human primary tumor tissues in the TCGA dataset ordered by similarity (n = 4). d Heatmap showing chromatin accessibility of the 7966 differentially accessible regions for the tIMEC, XD, PT, and MD samples (z-score of the log2-CPM is presented, cut-off log2fold change > 1 and FDR < 0.1) in the four biological replicates. e Density plots of H3K4me1, H3K27ac, and MED1 CUT&RUN signals in the four differential clusters identified by ATAC-seq in tIMEC, XD, and MD cells. Signal distribution is centered on the TSS peak and ranges from −5kB to 5kB. Greyscale signal intensity legend is valid for the three cell lines. The box plots in (c) indicate median values (middle lines), first and third quartiles (box edges), and 10th and 90th percentiles (error bars).

Fig. 3. Long-range interactions affect chromatin looping during tumor progression.

a H3K27ac HiChIP raw interaction maps of the MYC locus in tIMEC, XD, and MD cells (from top to bottom). 250 kb to 5 kb resolution (left to right). MYC locus (window = 80 Mb left, 3 Mb right). Numbers below the interaction maps correspond to maximum signal in the matrix. b Snapshot of the genomic locus surrounding the MYC locus (hg19, chr8:126,557,447-129,597,447), showing the ATAC-seq signal, insulation scores, H3K27ac interactions defined by HiChIP and identified TAD boundaries for tIMEC, XD, and MD. In addition, boundaries for TNBC cell lines (HCC1599 and MB157) and CTCF sites within these boundaries are shown. c Comparison of degree distributions (connection/region), according to type of connection: Promoter-Promoter (P-P), Enhancer-Enhancer (E-E), and Promoter-Enhancer (P-E). p-values have been computed, using the Kruskal–Wallis H-test. d Representation of the iCD of MEF2A Promoter (yellow) - Enhancer (gray) interactome in tIMEC (top), XD (middle), and MD (bottom). e Network average degree after MD-specific CREs removal (red line) versus network average degree after the removal of the same number of random regions (distribution) in tIMEC (left), XD (center), and MD (right) cells. f Kde plot, showing the relationship between the abundance of expressed genes (y-axis) and their level of expression (x-axis). Genes are clustered in 10 percentiles, according to their level of HiChIP connectivity. g Heatmap of differential RNA-seq transcripts z-scores. Genes are divided in three clusters, whose count is expressed on the right.

Fig. 4. TFs contributes to the activation of Metastatic-specific enhancer.

a Scheme of the ATAC-seq and RNA-seq data integration workflow using the IMAGE machine learning approach (created with BioRender.com). b Clustering of top MD-associated motifs with a motif activity MD > tIMEC and XD cells and defined as causal for gene regulation with an FDR < 0.01. Additional columns show z-score normalized expression of motif-associated TF, PCC between motif activity, and gene expression of associated TF and the associated TF family. Dot size for the TF family indicates enrichment in the unfiltered list by hypergeometric t-test. c Violin plots showing activity of MD-enriched motifs in tIMEC, XD, and MD cells, with an activity MD > tIMEC/XD. d Violin plots showing normalized expression values of motif-associated TFs showing a higher motif activity in MD compared to tIMEC or XD. e, f Distribution plot showing the rank of the TFs identified as drivers of MD enhancer activity among all TFs considered for the analysis (e) or expressed in the TCGA BC dataset (f). g Density plots showing the SOX9 CUT&RUN signal in the four differential ATAC-seq clusters identified. Signal distribution is centered on the TSS peak and ranges from −5kB to 5kB. Greyscale signal intensity legend is valid for the three cell lines. h Cumulative plots showing SOX9 CUT&RUN signal at peaks identified by ATAC-seq for the four defined clusters. i Representation of SOX9 iCD in MD cells. Vertices colors represent SOX9-linked enhancers (light blue), enhancers (gray), and promoters (yellow). The grayscale of edges is proportional to the number of connections between vertices. The violin plots in (c) and (d) indicate median values (middle lines), first and third quartiles (dashed lines). Statistical significance was determined by one-tailed unpaired student’s t-test.

Fig. 5. Enhancer redundancy specifies SOX9 transcriptional memory.

a Snapshot of the genomic locus surrounding the SOX9 (hg19, chr17:69400000-70700000), showing the ATAC-seq signal, the insulation scores, insulator boundaries for tIMEC, XD, and MD. In addition, CTCF sites for TNBC cell lines (HCC1599 and MB157), the 20 enhancers identified by ATAC-seq around the SOX9 locus and H3K27ac HiChIP enhancer-promoter (E-P) and enhancer-enhancer (E-E) interactions in MD. b Schematic representation of transcriptional memory in response to atRA. Cells are primed with a 2 h 1 µM atRA pulse, followed by a 36 h chase phase, and then treated again with 1 µM atRA or vehicle for 30 min or 1 h (created with BioRender.com). c Barplots of SOX9 relative expression in naïve and primed tIMEC, XD, and MD cells, untreated or stimulated with atRA for 30 or 60 min. Data are means of three independent biological replicates ± S.E.M. d Heatmap showing the changes in expression level (z-score of spike-in normalized FPKM values) identified by nascent RNA-seq in three biological replicates. e Uniform manifold approximation and projection (UMAP) embedding showing individual cells labeled according to their cell type (left), highlighting the enrichment of memory signature (right); number of analyzed cells: 22158. f Breadth distributions of H3K4me3 CUT&RUN signals in MD cells at the TSS of either all analyzed genes (dashed line) or memory genes cluster (blue line). On the top right, three genome browser snapshots are shown for H3K4me3 domains of different sizes from a memory (top), responsive (middle) or a non-responsive (bottom) gene. Window size for all three plots is 10 kb. g (Left) Loss of HiChIP regions associated with memory (blue) and responsive (yellow) promoters, upon filtering according to increasing thresholds of connectivity. (Right) Visualization of SOX9 clusters without connectivity filters (top), using a connectivity threshold equal to the 30th percentile of connectivity (middle) and using a connectivity threshold equal to the 40th percentile of connectivity (bottom). h Strength of chromatin annotations associated with memory (top) and responsive (bottom) clusters. The value of strength is obtained through the ratio between the fraction of HiChIP regions annotated with a particular chromatin state in non-filtered clusters and in clusters filtered with a threshold corresponding to the 90th percentile of connectivity. Statistical significance was determined by one-tailed unpaired student’s t-test.

Fig. 6. RA-mediated transcriptional memory relies on enhancer activity.

a Representative images of SOX9 nascent RNA detected by smRNA-FISH in naïve and primed MD cells, after 30 min treatment with vehicle or atRA. Nascent SOX9 RNA, red; DAPI, blue; scale bar = 5 µm. b Violin plots showing SOX9 nascent RNA foci area (pxl) in naïve and primed MD cells after treatment with vehicle or atRA. c Barplots of variation coefficient (CV) of SOX9 nascent RNA foci size. d Representative images of SOX9 nascent RNA detected by smRNA-FISH coupled with BRD4 immunostaining in naïve and primed MD cells after 30 min treatment with vehicle or atRA. DAPI, blue; BRD4, green; nascent SOX9 RNA, red; scale bar = 3 µm; smaller crop scale bar = 0.5 µm. e Violin plots showing percentage of colocalization between SOX9 nascent RNA foci and BRD4 immunostaining in naïve or primed MD cells, after treatment with vehicle or atRA. f Schematic representation of CRISPRi experimental setup obtained by combining the SunTag-dCAS9 with scFV-KRAB and BFP-sgRNAs targeting MD-enriched enhancers (#13, 14, 15 and 16), enhancers shared among tIMEC, XD and MD cells (#6 and 17), scramble (SCR), SOX9 promoter or a distal unrelated element (created with BioRender.com). g Violin plots showing the fold change of SOX9 nascent RNA foci area, between primed cells treated with atRA in respect to naïve condition, upon epigenome silencing of control (SCR), MD-enriched (b) or shared Enhancers (c). h Barplots of SOX9 relative expression levels in MD cells treated with vehicle or CDK8i at 10 nM or 100 nM. Comparison of vehicle and treatments is represented as fold change of expression of the target amplicon in primed cells after 30 min atRA treatment with respect to vehicle-treated MD cells. The barplots in (c) and (h) are means of three independent biological replicates ± S.E.M. The violin plots in (b), (e), and (g) indicate median values (middle lines), first and third quartiles (dashed lines) retrieved from three (b, e) or two (g) independent biological replicates. Statistical significance was determined by one-tailed unpaired student’s t-test (b, c, g, h) or one-tailed Mann–Whitney test (e).

Fig. 7. SOX9 drives metastatic-specific cell dormancy through RA-response.

a Representative images and quantification of p27 positive cells in MD expressing shGFP (Ctrl) or shSOX9 (KD SOX9) in steady state or after 72 h treatment with atRA. p27, white; DAPI, blue; scale bar = 20 µm. b Representative histograms of FACS analysis of Ctrl MD and KD SOX9 MD cells showing CellVue Maroon staining intensity (APC-A) 8 days (T1) or just after staining (T0), in the indicated conditions; barplot of the percentage of quiescent cell is shown. c Schematic representation of the experimental setup for NK cytotoxicity assays. (created with BioRender.com). d Average percentage of NK-mediated cytotoxicity of cycling and quiescent MD cells expressing either shGFP (Ctrl) or shSOX9 (KD SOX9), for different effector:target (E:T) ratios. e Barplots of relative gene expression levels of NK activating (green) and inhibitory (red) ligands in bulk and quiescent (p27 + ) MD cells. f Average percentage of NK-mediated cytotoxicity of quiescent MD cells for different effector:target (E:T) ratios, after vehicle or Navitoclax treatments. g Schematic representation of the experimental setup for metastatic colonization assay (created with BioRender.com). h DTCs and metastatic foci distribution within the lungs of nude mice injected with Ctrl or KD SOX9 MD cells. Top, tile scans of entire lung sections; scale bar: 1 mm. Bottom, representative images of metastatic foci (corresponding to i and ii from top image; scale bar = 1 mm) and quiescent DTCs (corresponding to iii-v from top image; scale bar = 50 µm; further zoom scale bar = 10 µm). mCherry: yellow; p27_Venus: green; DAPI: blue. i Barplots of scattered DTCs and metastatic foci in lungs of mice injected with Ctrl or KD SOX9, normalized to the lung area analyzed (n = 16; data combined from two experimental groups). j Barplots of the percentages of cycling or quiescent cells retrieved as scattered DTCs or lung metastases in nude mice (n = 16; data combined from two experimental groups). k Schematic representation of the experimental setup for metastatic colonization assay in humanized mouse models. (created with BioRender.com) l DTCs and metastatic foci distribution within the lungs of NSG mice injected with Ctrl or KD SOX9 MD cells and with or without NK-92 cells. Right, tile scans of entire lung sections; scale bar = 1 mm. Left, representative images of metastatic foci (corresponding to i and iii from top image; scale bar = 1 mm) and quiescent DTCs (corresponding to ii, iv, v, and v from top image; scale bar = 50 µm). Human mitochondria: yellow; p27_Venus: green; DAPI: blue. m Barplot of scattered DTCs and metastatic foci normalized to the lung lobe area analyzed in lungs of NSG mice injected with Ctrl or KD SOX9 MD cells, and with or without NK-92 cells, (n = 15; data combine from three experimental groups). n Barplots of the percentages of cycling or quiescent cells retrieved as scattered DTCs and metastatic foci in the lungs of NGS mice (n = 15; data combined from three independent experimental groups). The barplots in (a, b, e) are means of three independent biological replicates ± S.E.M. The mean values ± S.E.M in (d) and (f) were retrieved from three independent biological replicates. Statistical significance was determined by one-tailed unpaired student’s t-test.