The development of intestinal organoids from single adult intestinal stem cells in vitro recapitulates the regenerative capacity of the intestinal epithelium1,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.
This is a preview of subscription content
Subscription info for Chinese customers
We have a dedicated website for our Chinese customers. Please go to naturechina.com to subscribe to this journal.
Rent or Buy article
Get time limited or full article access on ReadCube.
All prices are NET prices.
All RNA-seq data generated here (bulk RNA-seq data (compound treatment at day 3), bulk RNA-seq data (compound treatment at day 0) and single-cell RNA-seq data) are available at the GEO under accession number GSE147136. The RXR antagonist used here (Cpd2170) is available through a material transfer agreement (MTA) with Novartis. The source data for following figures are available in the manuscript files: Figs. 1a, b, 2a, f, 3a–d, f, g and Extended data Figs. 1c–i, 2a–h, 3a–f, 4b–d, i, j, 5b, j–l, 6b, c, 7b, d–g, i, 8a–l, 9a, b, f, g, 10b, e, f, h. Gene sets used for visualization in Fig. 3e and Extended Data Fig. 9c–e, g are provided in Supplementary Table 2. Source data are provided with this paper.
Code used for image analysis was developed in the Liberali laboratory using MATLAB and Python 3. Segmentation for the image-based screen was performed using code developed in the Liberali laboratory in MATLAB, and is available at https://github.com/fmi-basel/glib-lukonin-et-al-2020. The code for organoid two-dimensional segmentation and feature extraction in other assays is available at https://github.com/fmi-basel/glib-nature2018-materials.
Serra, D. et al. Self-organization and symmetry breaking in intestinal organoid development. Nature 569, 66–72 (2019).
Grün, D. et al. Single-cell messenger RNA sequencing reveals rare intestinal cell types. Nature 525, 251–255 (2015).
Garabedian, E. M., Roberts, L. J., McNevin, M. S. & Gordon, J. I. Examining the role of Paneth cells in the small intestine by lineage ablation in transgenic mice. J. Biol. Chem. 272, 23729–23740 (1997).
Sato, T. et al. Paneth cells constitute the niche for Lgr5 stem cells in intestinal crypts. Nature 469, 415–418 (2011).
Marshman, E., Booth, C. & Potten, C. S. The intestinal epithelial stem cell. BioEssays 24, 91–98 (2002).
Bullen, T. F. et al. Characterization of epithelial cell shedding from human small intestine. Lab. Invest. 86, 1052–1063 (2006).
de Sousa e Melo, F. & de Sauvage, F. J. Cellular plasticity in intestinal homeostasis and disease. Cell Stem Cell 24, 54–64 (2019).
Tetteh, P. W. et al. Replacement of lost Lgr5-positive stem cells through plasticity of their enterocyte-lineage daughters. Cell Stem Cell 18, 203–213 (2016).
van Es, J. H. et al. Dll1+ secretory progenitor cells revert to stem cells upon crypt damage. Nat. Cell Biol. 14, 1099–1104 (2012).
Tian, H. et al. Opposing activities of Notch and Wnt signaling regulate intestinal stem cells and gut homeostasis. Cell Rep. 11, 33–42 (2015).
Joosten, S. P. J. et al. MET signaling mediates intestinal crypt-villus development, regeneration, and adenoma formation and is promoted by stem cell CD44 isoforms. Gastroenterology 153, 1040–1053 (2017).
Sato, T. et al. Single Lgr5 stem cells build crypt-villus structures in vitro without a mesenchymal niche. Nature 459, 262–265 (2009).
Chacón-Martínez, C. A., Koester, J. & Wickström, S. A. Signaling in the stem cell niche: regulating cell fate, function and plasticity. Development 145, dev165399 (2018).
Chang, J. et al. Proteomic changes during intestinal cell maturation in vivo. J. Proteomics 71, 530–546 (2008).
Carithers, L. J. et al. A novel approach to high-quality postmortem tissue procurement: the GTEx project. Biopreserv. Biobank. 13, 311–319 (2015).
Peeters, T. & Vantrappen, G. The Paneth cell: a source of intestinal lysozyme. Gut 16, 553–558 (1975).
VanDussen, K. L. et al. Notch signaling modulates proliferation and differentiation of intestinal crypt base columnar stem cells. Development 139, 488–497 (2012).
Yin, X. et al. Niche-independent high-purity cultures of Lgr5+ intestinal stem cells and their progeny. Nat. Methods 11, 106–112 (2014).
Levine, J. H. et al. Data-driven phenotypic dissection of AML reveals progenitor-like cells that correlate with prognosis. Cell 162, 184–197 (2015).
Liberali, P., Snijder, B. & Pelkmans, L. A hierarchical map of regulatory genetic interactions in membrane trafficking. Cell 157, 1473–1487 (2014).
Snijder, B., Liberali, P., Frechin, M., Stoeger, T. & Pelkmans, L. Predicting functional gene interactions with the hierarchical interaction score. Nat. Methods 10, 1089–1092 (2013).
Andreu, P. et al. A genetic study of the role of the Wnt/beta-catenin signalling in Paneth cell differentiation. Dev. Biol. 324, 288–296 (2008).
Pinto, D., Gregorieff, A., Begthel, H. & Clevers, H. Canonical Wnt signals are essential for homeostasis of the intestinal epithelium. Genes Dev. 17, 1709–1713 (2003).
Sakaki, J. et al. Synthesis and structure-activity relationship of novel RXR antagonists: orally active anti-diabetic and anti-obesity agents. Bioorg. Med. Chem. Lett. 17, 4804–4807 (2007).
Yui, S. et al. YAP/TAZ-dependent reprogramming of colonic epithelium links ECM remodeling to tissue regeneration. Cell Stem Cell 22, 35–49 (2018).
Gregorieff, A., Liu, Y., Inanlou, M. R., Khomchuk, Y. & Wrana, J. L. Yap-dependent reprogramming of Lgr5+ stem cells drives intestinal regeneration and cancer. Nature 526, 715–718 (2015).
Borel, P. et al. Processing of vitamin A and E in the human gastrointestinal tract. Am. J. Physiol. Gastrointest. Liver Physiol. 280, G95–G103 (2001).
Lampen, A., Meyer, S., Arnhold, T. & Nau, H. Metabolism of vitamin A and its active metabolite all-trans-retinoic acid in small intestinal enterocytes. J. Pharmacol. Exp. Ther. 295, 979–985 (2000).
Haber, A. L. et al. A single-cell survey of the small intestinal epithelium. Nature 551, 333–339 (2017).
Goverse, G. et al. Diet-derived short chain fatty acids stimulate intestinal epithelial cells to induce mucosal tolerogenic dendritic cells. J. Immunol. 198, 2172–2181 (2017).
Molotkov, A. & Duester, G. Genetic evidence that retinaldehyde dehydrogenase Raldh1 (Aldh1a1) functions downstream of alcohol dehydrogenase Adh1 in metabolism of retinol to retinoic acid. J. Biol. Chem. 278, 36085–36090 (2003).
Bowles, J. et al. ALDH1A1 provides a source of meiosis-inducing retinoic acid in mouse fetal ovaries. Nat. Commun. 7, 10845 (2016).
Ordóñez-Morán, P., Dafflon, C., Imajo, M., Nishida, E. & Huelsken, J. HOXA5 counteracts stem cell traits by inhibiting Wnt signaling in colorectal cancer. Cancer Cell 28, 815–829 (2015).
Nusse, Y. M. et al. Parasitic helminths induce fetal-like reversion in the intestinal stem cell niche. Nature 559, 109–113 (2018).
Gao, N., White, P. & Kaestner, K. H. Establishment of intestinal identity and epithelial-mesenchymal signaling by Cdx2. Dev. Cell 16, 588–599 (2009).
Spence, J. R. et al. Directed differentiation of human pluripotent stem cells into intestinal tissue in vitro. Nature 470, 105–109 (2011).
Kim, T. H. & Shivdasani, R. A. Genetic evidence that intestinal Notch functions vary regionally and operate through a common mechanism of Math1 repression. J. Biol. Chem. 286, 11427–11433 (2011).
Bhattacharya, N. et al. Normalizing microbiota-induced retinoic acid deficiency stimulates protective CD8+ T cell-mediated immunity in colorectal cancer. Immunity 45, 641–655 (2016).
Dotti, I. et al. Alterations in the epithelial stem cell compartment could contribute to permanent changes in the mucosa of patients with ulcerative colitis. Gut 66, 2069–2079 (2017).
Serafini, E. P., Kirk, A. P. & Chambers, T. J. Rate and pattern of epithelial cell proliferation in ulcerative colitis. Gut 22, 648–652 (1981).
Simmini, S. et al. Transformation of intestinal stem cells into gastric stem cells on loss of transcription factor Cdx2. Nat. Commun. 5, 5728 (2014).
Evans, R. M. & Mangelsdorf, D. J. Nuclear receptors, RXR, and the big bang. Cell 157, 255–266 (2014).
Cassani, B., Villablanca, E. J., De Calisto, J., Wang, S. & Mora, J. R. Vitamin A and immune regulation: role of retinoic acid in gut-associated dendritic cell education, immune protection and tolerance. Mol. Aspects Med. 33, 63–76 (2012).
Rangan, P. et al. Fasting-mimicking diet modulates microbiota and promotes intestinal regeneration to reduce inflammatory bowel disease pathology. Cell Rep. 26, 2704–2719 (2019).
Nepusz, T., Yu, H. & Paccanaro, A. Detecting overlapping protein complexes in protein-protein interaction networks. Nat. Methods 9, 471–472 (2012).
Bindea, G. et al. ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics 25, 1091–1093 (2009).
Barisic, D., Stadler, M. B., Iurlaro, M. & Schübeler, D. Mammalian ISWI and SWI/SNF selectively mediate binding of distinct transcription factors. Nature 569, 136–140 (2019).
We thank D. Vischi, S. Iftikhar and E. Tagliavini for IT support; L. Gelman for assistance and training; H. Kohler for sorting; S. Smallwood for sequencing; and L. Pelkmans, J. Betschinger, C. Tsiairis, S. Gasser and laboratory members for reading the manuscript. This work received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement number 758617) and from the Swiss National Foundation (POOP3_157531 to P.L.).
P.L. and I.L. are inventors on the patent application EP19182782.3, filed on 27 June 2019, with the title ‘Promoting tissue regeneration’, pertaining to the use of RXR antagonists as therapeutic agents in tissue regeneration. The patent applicant is the Friedrich Miescher Institute for Biomedical Research. Other authors declare no competing interests.
Peer review information Nature thanks Hans Clevers, Francesco Ioro, Bon-Kyoung Koo and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Extended data figures and tables
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; StdIntensity, 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.
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.
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. Senr, 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 high-compound 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.
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 co-annotated 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 CaCl2 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.
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. R2, 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 × 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.
a, Network representation of annotation enrichment analysis for a published RNA-seq dataset describing marker genes of intestinal epithelium cell types29. 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 log2-normalized transcripts per kilobase million (TPM) values) of Aldh1a1, Aldh1a2 and Aldh1a3 across cell types in the published single-cell RNA-seq (scRNA-seq) dataset29. 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 log2 normalized transcript counts) of Aldh1a1, 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 (log2 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.
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 (log2 fold change, log2FC). 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 × 10−3; ***P < 2 × 10−6. i, Left, expression of genes reported to have biased expression in the indicated gastrointestinal-tract tissues during unperturbed intestinal organoid development. log2-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 × 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, log2-transformed CPM values, normalized to mean expression over all samples.
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 log2CPM) 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 infection34. Expression is quantified by log2CPM, 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 (log2CPM) of genes differentially regulated in helminth infection34 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). log2CPM, 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. log2CPM, 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 × 10−23; DMSO versus atRA, 8.17 × 10−7; atRA versus RXRi, 5.84 × 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.
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 (log2 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. 29; n = 29,565 single cells. c–e, tSNE maps of our scRNA-seq experiment, colour coded by mean normalized expression (log2-transformed counts) of gene groups. c, Expression of cell-type signature genes from ref. 29. The numbers of genes per cell-type signature are shown in individual plots. d, Expression of genes that are differentially regulated following intestinal helminth infection34. 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 log2-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 log2-transformed fold change (log2FC) in the ‘early progenitor’ population compared with the ‘regenerative identity’ population. Every data point represents a gene. Centre, log2FC 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 (log2-transformed counts) of genes with log2FC 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 of z-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.
This file contains: Supplementary Methods and additional references and Supplementary Figure 1.
Description of the annotated compound library and screening hits This table includes (1) list of compounds comprising the screening library used in this study with corresponding p-values (see Supplementary Methods and Extended Data Fig.2), n=2242 independent conditions used for p-value calculation, (2) list of annotated primary targets of the compound library with gene-level p-values (Fisher’s exact test), n=895 target genes used for p-value calculation, (3) list of identified hits with respective annotated primary targets (related to Fig. 1b), (4) results of validation re-screen, n=49 compounds, (5) list of identified high confidence hits (related to Extended Data Fig. 3), (6) list of HIS-inferred functional interactions in tabular form, n=6158 interaction pairs, (7) list of input genes for calculation of HIS (related to Fig. 1b and Extended Data Fig. 4), n=207 genes.
Results of the compound treatment RNA sequencing experiments This table includes the normalized data as log2(cpm) mean centered and average of replica and differential expression in contrasts as log2 fold change (FC) values for the bulk RNA sequencing datasets with treatment at day 3 ((1) and (2), respectively, n=1884 genes) and at day 0 ((3) and (4) respectively, (3): n=4176 genes, (4): n=13527 genes), (5) summary of the scRNAseq experiment, n=29565 single cells (6) single cell-level analysis of heterogeneity in RXRi treated organoids (n=40 differentially regulated genes), (7) – (11) marker genes used in cross-correlation studies, n=1720, 652, 427, 16 and 593 genes respectively.
Antibodies and oligonucleotides used in this study This table includes (1) antibody stocks, dilutions and incubation times used for immunofluorescence staining in this study, (2) synthetic guide sequences used for generating lentiviral CRISPR-Cas9 constructs and corresponding oligonucleotide sequences.
About this article
Cite this article
Lukonin, I., Serra, D., Challet Meylan, L. et al. Phenotypic landscape of intestinal organoid regeneration. Nature 586, 275–280 (2020). https://doi.org/10.1038/s41586-020-2776-9
Experimental & Molecular Medicine (2021)
Nature Reviews Materials (2021)
Morphological screening of mesenchymal mammary tumor organoids to identify drugs that reverse epithelial-mesenchymal transition
Nature Communications (2021)
Cell Research (2021)
Experimental & Molecular Medicine (2021)