Phenotypic landscape of intestinal organoid regeneration

Abstract

The development of intestinal organoids from single adult intestinal stem cells in vitro recapitulates the regenerative capacity of the intestinal epithelium^1,2. Here we unravel the mechanisms that orchestrate both organoid formation and the regeneration of intestinal tissue, using an image-based screen to assay an annotated library of compounds. We generate multivariate feature profiles for hundreds of thousands of organoids to quantitatively describe their phenotypic landscape. We then use these phenotypic fingerprints to infer regulatory genetic interactions, establishing a new approach to the mapping of genetic interactions in an emergent system. This allows us to identify genes that regulate cell-fate transitions and maintain the balance between regeneration and homeostasis, unravelling previously unknown roles for several pathways, among them retinoic acid signalling. We then characterize a crucial role for retinoic acid nuclear receptors in controlling exit from the regenerative state and driving enterocyte differentiation. By combining quantitative imaging with RNA sequencing, we show the role of endogenous retinoic acid metabolism in initiating transcriptional programs that guide the cell-fate transitions of intestinal epithelium, and we identify an inhibitor of the retinoid X receptor that improves intestinal regeneration in vivo.

Extended Data Fig. 1. Reproducibility of the image-based screen.

a, An example image (MIP of a confocal z-stack) of an organoid, with annotated regions along the major axis, and showing antibody staining for lysozyme (Lyz) and aldolase B (Aldob). Lysozyme-positive Paneth cells are indicated by white arrowheads. Scale bar, 50 μm. b, Left, an example image of a stitched well-level overview. Scale bar, 100 μm. Middle, enlarged regions from the left panel. Scale bar, 100 μm. Right, segmentation masks, with object unique identifier numbers overlaid. a, b, Representative images of the DMSO control condition, repeated in n = 471 independent replicates with similar results. c, Top, reproducibility of detected organoid counts in replicas of the library plates. Heat map of Pearson’s correlation of organoid counts in screen replicates, with plate replicates grouped by parent library plates as indicated; n = 384 wells per plate. Bottom, heat map of standard deviation of organoid size (z-scored organoid area) for respective well positions. Passive and active controls are in wells outlined by the white box. n = 18 independent plates. d, Left, covariance matrix of the 34 extracted features, clustered by covariance (c-means clustering, ten clusters). Ch, channel; Stdlntensity, standard deviation of intensity. Right, heat map showing principal component loading of the extracted features; bar chart above, eigenvalues of first ten principal components. e, Top, changes in the accuracy of the naive Bayes predictor throughout the backwards-feature- elimination test. Bottom, final set of nine features (‘Algorithm-selected’) and list of features used as input for clustering with PhenoGraph. f, Distribution of mean values per well of the indicated features in active and passive control conditions. CHIR is the GSK3B inhibitor CHIR99021; DAPT is a γ-secretase inhibitor. Centre line, median; boxes, quartile range; individual data points overlaid. Two-sided t-test, with P-values shown in the plot. MeanInt, mean intensity. n = 471, 100 and 101 independent samples for the DMSO, CHIR99021 and DAPT conditions respectively. g, Z-scored organoid features used for clustering, shown for individual well positions of in-plate replicates of active and passive control conditions. n = 6, 6 and 20 wells per plate for DAPT, CHIR and DMSO respectively; shown are mean values over n = 18 assay plates. D, DAPI; Eccentr., eccentricity, Integr. Int., integrated intensity; Mass disp., mass displacement for the indicated immunofluorescence channels; Std Int., standard deviation of intensity. h, Mean z-scored values for the indicated features per well from the entire dataset, sorted in ascending order. Every data point is an average value per well, with conditions colour-coded as in the key. n = 5,250 independent samples measured in the screen. i, tSNE maps colour coded for z-scored values of indicated features (left to right: eccentricity, organoid area, mean intensity of aldolase B and mean intensity of lysozyme antibody staining) and PhenoGraph classification. Every data point represents an individual organoid, with sampling stratified by PhenoGraph assignment. The sample size is 5,000 individual organoids.

Extended Data Fig. 2. Robust phenotypic signatures.

a, Top, phenotypic clusters and phenotypic groups (see Fig. 1b) identified using the multivariate feature array, and heat map of mean z-scored values per class for the indicated feature array. Aldob, aldolase B; Ar, area; D, DAPI; Ecc, eccentricity; II, integrated intensity; Lyz, lysozyme; MD, mass displacement; MI, mean intensity; SI, standard deviation of intensity; Sh, shape. n = 40,624, 34,265, 20,134, 31,687, 2,549, 702, 23,643, 48,992, 30,791, 21,342, 23,610, 32,773, 19,172, 41,931 and 29,939 organoids for classes 1-15 respectively. Bottom, phenotypic class assignment for organoids in control conditions; outer circle colour, phenotypic group; inner circle colour, assigned phenotypic cluster; n = 7367, 7982 and 33395 individual organoids for DAPT, CHIR99021 (CHIR) and DMSO conditions respectively. b, Potential loss of phenotypic resolution owing to phenotype grouping. Three phenotypic clusters are assigned to the ‘regenerative’ phenotype, but the exact 15-element phenotypic signature differs depending on whether PORCN or GSK3β is inhibited. Mean of n = 100, 4, 2 and 2 independent replicates for CHIR99021, Cpd2170, Cpd420 and Cpd1772 respectively. c, Top, relative z-score transformed abundances of 15 phenotypic clusters, shown as means ± s.d. Bottom, heat maps showing individual and mean phenotypic signatures for 50 replicates of the indicated control conditions. n = 50 sampled independent replicates for all 3 conditions. d, Left, randomization trial for reshuffling of phenotypic class labels in the entire dataset. Right, scatter plots and histograms showing the distribution of replica correlation and maximum class variance in permuted and non-permuted data for the selected and filtered conditions. Red dashed lines show the original selection stringency. e, Distribution of calculated P-values in the indicated groups of conditions. Top, scatter plot and kernel density estimation; bottom, box plots. Boxes, quartile range; whiskers, value interval with excluded outliers; solid lines, median values. n = 2,242 unique conditions corresponding to 5,204 individual wells. n = 279, 15 and 1,948 conditions for ‘hit’, ‘hit (one replica)’ and ‘not selected’ groups respectively. f, Top left, mean z-scored values per identified phenotypic class for the indicated feature array. Bottom left, heat map representation of the phenotypic signature for the indicated compounds, grouped by the most abundant phenotype in the original screen. Top right, tSNE map colour coded for the highest enriched phenotypic clusters. Every data point represents a compound. tSNEs were calculated from phenotypic signatures. Bottom, calculated reproducibility values per compound and numbers of replicates used for analysis of the indicated compounds. The replica count is the number of non-sparse (n > 10 organoids) independent wells recorded per condition. Assay performed in n = 8 independent replica from 2 independent organoid cultures. g, Flowchart depicting filtering of the screened conditions. h, tSNE maps colour coded for the highest enriched phenotypic clusters as indicated; every data point represents either a compound condition (top; n = 301 compound conditions) or a unique target gene-modality pair (bottom; n = 207 target genes); tSNEs calculated from 15-dimensional phenotypic signatures.

Extended Data Fig. 3. Identification of hit compounds and target genes.

a, Top, compound abundance in the library and hit list. Bottom, distribution of enrichment scores for targets of the screen hits. Every data point represents a unique target-gene–modality pair. b, Heat map showing phenotypic signatures for genes with a compound coverage of at least three compounds per gene. Heat map rows marked by arrowheads show compounds used as representative conditions per target. S_enr, target enrichment score (see a). c, Left, contingency table for two-sided Fischer’s exact test; right, bar plot and box plot of P-value distributions in the indicated compound condition groups. High conf. hits, high-confidence hits. In the box plot, boxes show quartile ranges; whiskers show value intervals with excluded outliers; centre lines show median values. n = 699, 167 and 29 unique target genes for non-hits, other hits and high-confidence hits respectively. d, Enrichment and false discovery rate for functional annotation enrichments reported by the STRING database tool (https://string-db.org/) for the list of 29 target genes (e). e, Top left, Venn diagram showing the distribution of target genes by compound coverage. Main panel, network of functional interactions for the set of genes with highcompound coverage and hitlist specific enrichment. Edges show interactions known in the STRING database; thickness and line type represent the STRING combined score (see network key); nodes show unique target genes; node outlines show phenotypic clusters with the highest enrichment; circles indicate groups of similar functional annotations. f, Enrichment of KEGG terms (percentage of genes annotated with KEGG terms) in conditions from the hit list where the highest enriched phenotype is either ‘regenerative’ or ‘enterocyst’. Values are sorted by corresponding P-values (two-sided test based on hypergeometric distribution with Bonferroni correction). n = 59 and 61 genes for the enterocyst and regenerative phenotypes, respectively.

Extended Data Fig. 4. Validation of HIS predictive power and screen hits.

a, Principle of the hierarchical interaction score (HIS). HIS is calculated pairwise for a set of genes described by multivariate phenotypic signatures (left), to generate a matrix of graded interactions (middle) in which every interaction contains information on which readout contributed the most. In the resulting HIS-based network (right), the strength of the interactions (the thickness of the edges) reflects the similarity of the observed enrichment or depletion, and directionality is from the gene with the broader set of effects to the gene with the narrower set of effects. b, Heat maps of the HIS pairs for n = 207 target genes; interactions with non-zero HIS values are shown colour-coded for the phenotypic class that contributes the most to HIS inference (left) or for the HIS value (right). c, Top left, number of retained nodes; bottom left, number of STRING-reported and predicted edges (see Supplementary Methods); centre, fold enrichment of reported STRING interactions over prediction; right, percentage of co-annotated gene pairs with increasing HIS threshold cutoff (top right, KEGG co-annotations; bottom right, GO co-annotations (levels 5–15)). Red dashed line, HIS threshold (0.175); arrowheads, inflection points. Top and bottom right plots: blue line, polynomial fit. d, Percentage of coannotated gene pairs in indicated sets of genes in measured and randomized data (top, KEGG co-annotations; bottom, GO co-annotations (levels 5–15)). The corresponding number of genes is drawn either from the genome or the hitlist pool, as indicated. e, HIS-based interaction network for β-catenin. Node colour shows the highest enriched phenotypic class; edge colours show the phenotypic class contributing the most to HIS value inferences (see key); edge arrow shapes show HIS-inferred directionality. All edges directed to/from β-catenin are shown. f, Inducible NICD–GFP and HES1–GFP knock-in organoids with and without doxycycline (dox) treatment from day 0, showing loss of Paneth cells in NICD-overexpressing organoids. Nuclei are stained with DAPI (blue); lysozyme is stained with antibody (red). Scale bar, 50 μm. g, Day 5 puromycin-selected lentiCRISPRv2-mediated Casr knockout (CasrKO) organoids, with corresponding infection and selection controls (pLenti–EF1–EGFP organoids). Organoids were cultured from single cells in the indicated conditions. The upper six images show the perturbed crypt phenotype and lack of phenotypic response to CaCl₂ treatment in CasrKO organoids. Bottom, lentivirally introduced EGFP expression after puromycin selection. Scale bars, 30 μm. h, Top two rows, day 3 and day 5 puromycin-selected control (pLenti–EF1–EGFP) organoids; bottom row, lentiviral Akt1-overexpressing (Akt1 OE) organoids cultured from single cells, showing perturbed phenotypes and lack of enterocyte differentiation. Scale bars, 30 μm. i, Top, organoids from DMSO control and RXR antagonist conditions analysed in the screen. Scale bar, 50 μm. Bottom, corresponding phenotypic signatures. n = 2, 2 and 1 non-sparse (n > 10 organoids) independent replica for Cpd2170, Cpd2390 and Cpd2173, respectively. j, Abundance of phenotypic cluster 13 in hit conditions with a reproducibility of more than 0.5. Compounds targeting RXR and RAR are depicted with coloured dots. Histograms show distributions of respective values. n = 301 independent conditions. Microscopy images shown in f-i are MIPs of confocal z-stacks. f-h, Assays performed in n = 12 independent replica from 2 independent organoid cultures with similar results.

Extended Data Fig. 5. Profiling the effects of treatment with an RXR antagonist.

a, Images of organoids cultured in the presence of an RXR agonist (RXRa) in the indicated treatment conditions. Scale bar, 50 μm. Assay performed in n = 12 (DMSO) or n = 6 (other conditions) independent replica from 2 independent organoid cultures. b, Top, distribution of mean aldolase B staining intensity in control organoids and those treated at day 4. Bottom, abundance of the ‘enterocyst’ phenotype in the indicated conditions. Error bars show standard deviations between n = 6 (DMSO) or n = 3 (other conditions) independent replica. Assays were performed with the same number of replicas from two independent organoid cultures with similar results. c, Images of organoids fixed at indicated time points, showing antibody staining for lysozyme. Scale bar, 20 μm. d, Top two rows, images of organoids from Fucci2 mice, expressing a dual-colour cell-cycle reporter (mVenus-hGem(1/110), G2/M-phase reporter; mCherry–hCdt1(30/120), G1-phase reporter) cultured with or without RXR antagonist (RXRi), fixed at day 4. Top row, expression of mVenus–hGem(1/110) (FUCCI (G2/M)), mCherry–hCdt1(30/120)(FUCCI (G1)); bottom row, antibody staining for YAP1. Scale bar, 100 μm. Bottom row, enlargements of those areas outlined in the middle row. Scale bar, 20 μm. e, Images of organoids following 10 μM atRA treatment at day 3. Top row, scale bar, 50 μm. Bottom row, scale bar, 10 μm. f, Organoids fixed at day 5, 48 h after switch to ENR medium. Scale bar, 100 μm. g, Top, distribution of mean SOX9, and bottom, ki67 staining intensity in control and treated organoids. Organoids were fixed at day 5, 48 h after switching to ENR medium. Violin ranges extended by one standard deviation. Two-sided t-test; P-values shown in the plot. h, Images of wild-type organoids, treated as indicated at day 3, and fixed at day 5. Scale bar, 100 μm. i, Images of organoids expressing RARE–dsGFP at day 4. dsGFP, destabilized GFP; scale bar, 10 μm. Assay performed in n = 12 independent replicates from n = 3 independent organoid cultures with similar results. j, Histological images of mouse small intestinal epithelium from RARE-LacZ mouse. Regions along the crypt–villus axis are indicated. A single crypt flanked by villus regions, outlined in the left panel, is shown in blue at the right. Scale bar, 50 μm. Experiment performed in n = 4 independent histological samples from n = 2 independent mice with similar results. k, Left, RARE–dsGFP organoids in the indicated treatment conditions. Scale bar, 60 μm. Right, DAPI-normalized sum intensity (norm. int.) of RARE–dsGFP reporter and aldolase B in individual organoids. Data points represent individual organoids, colour-coded according to treatment as indicated. R ², Pearson’s correlation coefficient. l, Left, puromycin-selected ALDH1A2- knockout organoids. Scale bar, 50 μm. Right, ALDH1A1 staining intensity in control and ALDH1A1-knockout organoids. n = 36 and 17 organoids for control and ALDH1A1-knockout conditions respectively. ***P = 5 x 10^-18, two-tailed t-test. m, Images of organoids treated as indicated from day 0. ALDH1A1i, 5 μM ALDH1A1 inhibitor; atRA, 10 μM all-trans retinoic acid. Scale bar, 50 μm. Assay performed in n = 8 independent replica. n, Images of organoids treated with RXR antagonist from day 0 in medium with or without vitamin A. Organoids were fixed at day 5, 48 h after switch to ENR medium. Scale bar, 50 μm. o, Wild-type organoids cultured in medium with or without vitamin A. Scale bar, 100 μm. p, Organoids cultured in medium with or without vitamin A and fixed at the indicated time points. Scale bar, 20 μm. a, d–f, h, k–o, Microscopy images are composite MIPs of confocal z-stacks. c, i, p, Microscopy images are composite individual z-planes from the organoid middle plane. b, g, k, n indicates the number of individual organoids used as data points. b, l, Boxes show quartile ranges; whiskers show value intervals with excluded outliers; middle lines or white dots show median values. c–f, h, Assays performed in at least four independent replicates from at least two independent organoid cultures with similar results. a, c–f, h, i, k–p, Organoids were cultured in the indicated conditions from single cells.

Extended Data Fig. 6. Retinol metabolism in enterocytes.

a, Network representation of annotation enrichment analysis for a published RNA-seq dataset describing marker genes of intestinal epithelium cell types. Small nodes represent genes; large nodes represent functional annotations; edge thickness show gene-annotation assignments (grey) or term–term relations (black). Colours indicates cell-type specificities of annotations (see key). b, Box plots showing expression (in log₂-normalized transcripts per kilobase million (TPM) values) of Aldhlal, Aldh1a2 and Aldh1a3 across cell types in the published single-cell RNA-seq (scRNA-seq) dataset. Boxes show quartile ranges; whiskers represent value intervals with excluded outliers; middle lines show median values. Each data point represents a single cell. EC, enterocyte; TA, transient-amplifying. n = 1,522 single cells. c, Expression (in log₂ normalized transcript counts) of Aldhlal, Aldh1a2 and Aldh1a3 across indicated treatment conditions in scRNA-seq (Fig. 3 and Extended Data Fig. 9). Sampling of n = 3,000 single cells. d, tSNE maps of the scRNA-seq experiment. Each data point represents a single cell, colour-coded by expression (log₂ normalized transcript counts) of indicated genes in samples of the scRNA-seq experiment (Fig. 3). Ranges of the colour mapping are indicated in each plot. Sampling of n = 5,000 single cells. b–d, Each data point represents a single cell. atRA, 10 μM all-trans retinoic acid; DMSO, DMSO control, RXRi, 5 μM RXR antagonist; w/o VitA, vitamin-A-depleted medium.

Extended Data Fig. 7. Transcriptional signature of organoids treated with RXR antagonist.

a, Motif enrichments in the 1,000-base-pair regions that are centred on promoters of genes that are more highly expressed at a given day compared with the average from the whole period of organoid development. b, Enrichment heat map of transcription-factor-binding motifs at the transition from day 3 to day 4. Genes are binned by expression pattern. c, Network of functional annotation (KEGG) terms specifically enriched in genes containing RXR- and RAR-containing transcription-factor-binding motifs. Nodes show annotation terms; edges represent term–term interactions; node colours show annotation groups (see key). d, Principal component (PC) analysis of RNA-seq samples from the indicated treatment conditions and time points. The scatter plot shows the first two principal components. n = 4 (2 independent organoid cultures each from 2 independent mice) for all individual conditions. e, Pearson’s correlation values between samples of the time-course RNA-seq experiment and compound-stimulation RNA-seq experiment. n = 12,025 genes detected in both datasets. f, Left, expression of selected genes in untreated organoids at day 3 and in the indicated treatment conditions 24 h after treatment. Every data point represents a gene, colour-coded according to assigned cluster label (k-means clustering). Stratified sampling of 1,884 genes; n = 50 genes per cluster. FDR, false discovery rate. Top right, abundance of genes containing RXR/RAR-binding motifs in assigned clusters; bottom right, mean expression of genes from respective clusters in the RXR-antagonist condition compared with the DMSO control (log₂ fold change, log₂FC). Data shown as means ± s.d. n = 152, 102, 196, 361, 396, 109, 113, 290 and 165 for expression-profile clusters 1–9 respectively. g, tSNE maps colour-coded according to assigned cluster labels (k-means clustering). Left, tSNEs calculated from expression values measured in the bulk RNA-seq experiment (e). Right, genes showing treatment-specific upregulation and downregulation in the orthogonal compound condition. Every data point represents a gene. h, Annotation enrichment analysis (GO terms) for gene categories described in g. Asterisks depict statistical significance determined with two-tailed hypergeometric test with Bonferroni correction. *P < 0.05; **P < 2 x 10^-3; ***P < 2 x 10^-6. i, Left, expression of genes reported to have biased expression in the indicated gastrointestinal-tract tissues during unperturbed intestinal organoid development. log₂-transformed count per million (CPM) values are normalized to mean expression over all samples. Solid lines show mean expression per time point; opaque intervals show standard deviation. n = 83, 20, 63, 156, 259, 45 and 31 genes for colon (sigmoid), colon (transverse), oesophagus (junction), oesophagus (mucosa), not specific, small intestine and stomach respectively. Middle, expression of genes specific to the indicated gastrointestinal tissues in RXRi-treated organoids and DMSO controls after 24 h of treatment. Right, expression of small-intestine-specific genes in DMSO- and RXRi-treated organoids. Boxes show quartile ranges; whiskers denote value intervals with excluded outliers; middle lines show median values. ***P < 2 x 10^-19 (two-tailed t-test). d–i, atRA, 10 μM all-trans retinoic acid; 9cisRA, 10 μM 9-cis-retinoic acid; DMSO, DMSO control; RXRi, 5 μM RXR antagonist. g, i, log₂-transformed CPM values, normalized to mean expression over all samples.

Extended Data Fig. 8. Organoids treated with an RXR antagonist lose intestinal identity.

a, Principal component analysis of RNA-seq samples (treated at day 0) from different treatment conditions and time points. The scatter plot shows the first two principal components. n = 4 (2 technical replicates of 2 biological replicates) for all individual conditions. The bar chart above depicts the percentage of variance explained by the first eight principal components. b, Pearson’s correlation values between samples of the time- course RNA-seq experiment. n = 13,502 genes detected in all samples. c, Normalized gene expression (in log₂CPM) of the indicated genes in samples of the time-course RNA-seq experiment (Fig. 3b). Enteroc., enterocyte specific; G.C., goblet cell specific. d, Expression of genes from the indicated categories in treatment conditions at days 1 and 5 of organoid development. Data points represent genes. n = 14 and 16 genes for intestinal cell type markers and fetal- like and YAP target genes respectively (gene names are in b). e, Normalized gene expression in contrast (RXRi-treated over DMSO control, at day 4) for samples treated at days 0 and 3. Left and middle, colours show cell-type markers (left) or indicated gene categories (middle); dashed lines are at fourfold enrichment or depletion; opaque areas cover genes not significantly perturbed in either contrast; data points represent genes; colours represent indicated gene categories; marker types highlight the genes of interest. Right, only those genes that were upregulated on transition from day 0 to day 1 in the DMSO control are shown. Red dashed lines are at fourfold enrichment. Marker types highlight the genes of interest. f, Expression of genes grouped by response to intestinal helminth infection. Expression is quantified by log₂CPM, normalized to the mean expression over all samples in the respective datasets. Solid lines show mean values per time points; opaque intervals show standard deviations; dashed lines show fourfold enrichment or depletion. g, Normalized expression (log₂CPM) of genes differentially regulated in helminth infection³⁴ at day 4 (top) and day 5 (bottom). Each data point represents a gene; n = 141 and 286 genes for infection-upregulated and -downregulated genes respectively. h, Expression of tissue-specific genes over time points and treatments (DMSO, atRA and RXRi). log₂CPM, normalized to mean expression over all samples; dashed lines show fourfold enrichment or depletion. n = 63, 156, 259, 45 and 31 genes for oesophagus (junction), oesophagus (mucosa), not specific, small intestine and stomach respectively. i, Top, expression of tissue-specific genes (n = 156, 259, 45 and 31 genes for oesophagus (mucosa), not specific, small intestine and stomach respectively), and bottom, expression of genes regulated by Cdx2 (n = 340, 170 and 11,378 genes for Cdx2-upregulated, Cdx2-downregulated and not affected respectively) in unperturbed organoid development. log₂CPM, normalized to mean expression over all samples. j, Left, fold changes for the indicated genes in the indicated contrasts over the time-course RNA-seq experiment (Fig. 3b). Right, normalized gene expression in the indicated contrasts for samples treated at days 0 and 3; dashed lines show fourfold enrichment or depletion; the opaque area covers genes not significantly perturbed in either contrasts; colours show the indicated gene categories; every data point represents a gene, n = 340, 170 and 11,378 genes for Cdx2-upregulated, Cdx2- downregulated and not affected respectively. k, Expression of genes regulated by Cdx2 in organoids over the time course of the experiment. Solid lines show mean values per time point; opaque intervals show standard deviations; dashed lines show fourfold enrichment or depletion. l, Left panels, wild-type organoids at day 3 cultured from single cells; images are composite MIPs; scale bar, 10 μm. Right, intensity of Cdx2 staining; boxes show quartile ranges; whiskers show value intervals with excluded outliers; solid lines show median values. Asterisks indicate statistical significance; two-sided t-test; P-values are as follows: DMSO versus RXRi, 2.04 x 10^-23; DMSO versus atRA, 8.17 x 10^-7; atRA versus RXRi, 5.84 x 10^-24; n, number of organoids analysed in each condition. m, Cdx2-EGFP organoids at day 4, cultured from single cells. Images are composite MIPs. Scale bars, 50 μm. The right two columns show YAP1 and Cdx2 staining from respective composite images. Assay performed in n = 16 independent replicates from n = 2 independent organoid cultures with similar results. f, h, i, k, Solid lines show mean values per time points; opaque intervals show standar deviations. atRA, 10 μM all-trans retinoic acid; DMSO, DMSO control; RXRi, 5 μM RXR antagonist.

Extended Data Fig. 9. scRNA-seq analysis of phenotypes produced by RXRi treatment.

a, Left, tSNE maps from scRNA-seq. Each data point represents a single cell, colour-coded by sample (see key). n = 4,787, 3,570, 5,365, 4,963, 5,576 and 5,304 single cells for day 4 DMSO, RXRi, atRA and vitamin A-free, and day 1 DMSO control samples 1 and 2 respectively. Right, expression (log ₂normalized transcript counts) of indicated genes in samples of the scRNA-seq experiment (Fig. 3); sampling of n = 3,000 cells. b, Cell-type composition of scRNA-seq samples. Confidence cutoffs for cell-type assignment were set individually for each cell type on the basis of the mean expression levels of marker genes from ref. ; n = 29,565 single cells. c-e, tSNE maps of our scRNA- seq experiment, colour coded by mean normalized expression (log₂-transformed counts) of gene groups. c, Expression of cell-type signature genes from ref. . The numbers of genes per cell-type signature are shown in individual plots. d, Expression of genes that are differentially regulated following intestinal helminth infection. e, Expression of genes with biased expression in tissues (GTEX database; https://gtexportal.org/home/). f, Left, partial tSNE maps (top) and density map (bottom) of single cells from the RXR-antagonist-treated sample of the scRNA-seq experiment. In the density map, the x- and y-axes show the mean expression (mean of log₂-transformed counts) of genes from the indicated categories. Colours show cell counts in respective bins. Red dashed lines show cutoffs used to assign populations. Right, tSNE maps of the scRNA- seq experiment, showing assigned identities in the RXR-antagonist-treated sample. g, Left, average expression of genes in populations of cells as assigned in f. Colours depict the log₂-transformed fold change (log₂FC) in the ‘early progenitor’ population compared with the ‘regenerative identity’ population. Every data point represents a gene. Centre, log₂FC values of the indicated genes. Green, absorptive-lineage-specific genes; blue, secretory-lineage- specific genes. Right, tSNE maps of the scRNA-seq experiment colour-coded by mean normalized expression (log₂-transformed counts) of genes with log₂FC values of greater than 0.5 between ‘early progenitors’ and ‘regenerative’ populations. a, b, atRA, 10 μM all-trans retinoic acid; DMSO, DMSO control; RXRi, 5 μM RXR antagonist; w/o VitA, vitamin-A-depleted medium. a–f, In tSNE maps each data point represents a single cell; sampling of 5,000 single cells. c–e, g, Ranges of colour mapping and numbers of genes per category (n) are indicated in each plot.

Extended Data Fig. 10. Treatment with an RXR antagonist improves intestinal regeneration in an irradiation-induced mouse model of colitis.

a, irradiation-induced colitis model. PK, pharmacokinetic analysis. b, Change in body weight in the indicated treatment cohorts over the time frame of the study. Data shown are mean ± s.d values for cohorts of n = 6 independent mice per cohort. Non-irr., non-irradiated vehicle-treated cohort. c, Images of spleen extracted from mice from the indicated treatment cohorts. Scale bar, 1 cm. d, Small-intestine’ Swiss rolls’ from mice belonging to the indicated treatment cohorts and time points. Scale bar, 5 mm. e, Top, haemoccult score of stool samples from the indicated treatment cohorts; middle lines show mean; error bars show ± s.d.; n = 6 independent animals per cohort. Bottom, representative images of the haemoccult samples corresponding to the quantitative haemoccult scores. f, Serum levels of the RXRi compound at 7 h post-treatment at the indicated time points in RXRi-treated mice. Middle lines show means; error bars show ± s.d.; n = 6 independent animals per cohort. g, Top, small-intestine Swiss roll, showing the ileum, duodenum and jejunum. Scale bar, 5 mm. Middle, jejunum part of the intestine. Scale bar, 250 μm. Bottom, enlarged view of villi units; individual are villi indicated with arrowheads; the shaded area indicates measured villi length. Scale bar, 250 μm. h, Top, abundance of goblet cells, measured by staining with alcian blue, in histological samples at day 6 in mice from the indicated groups. Bottom, quantified percentage of alcian blue (AC)-positive cells in the entire intestine of mice from the indicated treatment cohorts. Boxes show, quartile ranges; whiskers show value intervals; solid lines show median values; individual data points are overlaid. n = 6 independent animals per cohort. i, Histological profiling of the abundance of Ki67+ cells in the intestinal epithelium at days 2, 4 and 6 in the indicated treatment cohorts. j–l, Histological profiling of the intestinal epithelium at days 2 (j), 4 (k) and 6 (l) extracted from mice from the indicated treatment cohorts. Left panels, haematoxylin and antibody staining for OLFM4; scale bars, 50 μm. Enlarged regions: upper insets, scale bars, 25 μm; lower insets, scale bars, 12 μm. Middle panels, nuclear fast red (NFR) and alcian blue staining; scale bars, 50 μm. Right panels, confocal imaging of histological sections (composite MIPs ofz-stacks); scale bars, 50 μm. Enlarged regions: upper insets, scale bars, 25 μm; lower insets, scale bars, 12 μm. In all panels, experiments were carried out in n = 6 independent mice per treatment cohort with similar results. c, d, Representative images from respective treatment cohort animals. g–l, Representative images from n = 24 histological samples from n = 6 independent mice per treatment cohort. b–h, Day 6 (terminal time point) mouse study repeated twice with n = 6 mice per cohort with similar results.

Fig. 1. Phenotypic landscape and map of functional genetic interactions in intestinal organoids.

a, Screening workflow: organoid-level multivariate feature array, clustering by similarity and generation of phenotypic fingerprints for each condition. D, DAPI; Aldob, aldolase B; Lyz, lysozyme; mass disp., mass displacement for the indicated immunofluorescence channels; tSNE, t-distributed stochastic neighbour embedding. b, Left, main panel: identified phenotypic clusters, representative images and schematic representation of the phenotypes. Right, network of HIS interactions (for HIS values of at least 0.175) between 101 genes. For the inner network, node colour and edge colour are explained in the key on the left. Highlighted nodes show off-diagonal interactors. The outer networks show detected subnetworks. Edge thicknesses show HIS values; edges with the highest HIS value per node are shown. Node outlines show compound modality (inhibitor or activator). KEGG enrichments are shown per phenotype (inner circle) and subnetwork (outer circle). Infl. med. reg., inflammation-mediated regulation. a, b, Images are maximum intensity projections (MIPs) of confocal z-stacks, showing nuclear staining (DAPI) and antibody staining for lysozyme and aldolase B.

Fig. 2. Retinol metabolism and RXR activity in enterocyte differentiation.

a, Top, images of organoids treated at day 0 as indicated (DMSO, dimethyl sulfoxide; control; RXRi, treatment with RXR antagonist Cpd2170). Bottom, quantified aldolase B intensity. b, Colour composite images show antibody staining for DLL1 and YAP1 under different treatments at the indicated time points. Scale bars, 20 μm. D, DAPI. Greyscale images show YAP1 localization in the highlighted regions. Scale bars, 10 μm. c, Upper four images show antibody staining for lysozyme, aldolase B and YAP1 under different treatments at day 3. Scale bar, 50 μm. Bottom, YAP1 localization in highlighted regions. Scale bar, 10 μm. d, Retinol metabolism. Cyp26, cytochrome P450; RE, retinyl ester. e, Effect of knocking out Aldhlal on organoid development. Top, images show antibody staining for ALDH1A1; bottom, antibody staining for lysozyme and aldolase B. ALDH1A1i, ALDH1A1 inhibition with A37; pLenti-EF1-EGFP, infection and selection control. f, Left, day 4 organoids under different treatments (ALDH1A1i, day 0; atRA, day 3), in medium with (control) and without vitamin A. Right, quantified aldolase B intensity. Asterisks show statistical significance for the same treatment with/without vitamin A; two-sided t-test; NS, not significant; P values shown in the plot. a-c, e, f, Representative images of organoids (n = 6 independent replicates from n = 3 independent experiments), cultured from single cells. b, Day 2 plus 9 h: organoid middle plane, day 3 and 4: MIP of confocal z-stacks. Nuclei are stained with DAPI; scale bar, 50 μm. a, f, Violin ranges extended by standard deviations; box plots show quartile ranges; whiskers show value intervals; white dots show median values; n values refer to the number of individual organoids in conditions.

Fig. 3. RXRi treatment prevents exit from regenerative state and improves regeneration in vivo.

a, Temporal expression of genes containing RAR or RXR motifs, normalized to expression at day 0. Solid lines show mean values; opaque interval show standard deviation; n values refer to the number of genes per group. b, Pearson’s correlations between gene-expression levels in time-course and compound experiments. c, Differential gene expression in organoids relative to DMSO at indicated time points. GC, goblet cell. d, Fold change relative to DMSO (day 4) for specific gene categories and cell-type markers at days 0 and 3. Data points represent individual genes; horizontal lines show fourfold enrichment or depletion; n, number of genes. e, tSNE maps generated from single-cell RNA sequencing data, colour-coded by sample (left) and by mean normalized expression of gene groups (right). Data points show single cells; sampling of 5,000; n, number of genes; colour map ranges are shown in plots. f, Change in body weight (AUC, area under curve) in treatment cohorts; n = 6 mice per cohort. g, Quantification of villus length; n, number of villi; measured in n = 3 independent samples per condition. In the images, tissues are stained with haematoxylin and eosin. a, c, d, log₂-transformed fold change (log₂(FC)). b–d, n = 4 independent cultures from 2 independent mice. d, f, g, Boxes show quartile range; whiskers show value intervals; solid lines show median values. Two-sided t-test; asterisks indicate statistical significance; individual P values shown in plots.