Harnessing single-cell genomics to improve the physiological fidelity of organoid-derived cell types

Harnessing single-cell genomics to improve the physiological fidelity of organoid-derived cell types Background: Single-cell genomic methods now provide unprecedented resolution for characterizing the component cell types and states of tissues such as the epithelial subsets of the gastrointestinal tract. Nevertheless, functional studies of these subsets at scale require faithful in vitro models of identified in vivo biology. While intestinal organoids have been invaluable in providing mechanistic insights in vitro, the extent to which organoid-derived cell types recapitulate their in vivo counterparts remains formally untested, with no systematic approach for improving model fidelity. Results: Here, we present a generally applicable framework that utilizes massively parallel single-cell RNA-seq to compare cell types and states found in vivo to those of in vitro models such as organoids. Furthermore, we leverage identified discrepancies to improve model fidelity. Using the Paneth cell (PC), which supports the stem cell niche and produces the largest diversity of antimicrobials in the small intestine, as an exemplar, we uncover fundamental gene expression differences in lineage-defining genes between in vivo PCs and those of the current in vitro organoid model. With this information, we nominate a molecular intervention to rationally improve the physiological fidelity of our in vitro PCs. We then perform transcriptomic, cytometric, morphologic and proteomic characterization, and demonstrate functional (antimicrobial activity, niche support) improvements in PC physiology. Conclusions: Our systematic approach provides a simple workflow for identifying the limitations of in vitro models and enhancing their physiological fidelity. Using adult stem cell-derived PCs within intestinal organoids as a model system, we successfully benchmark organoid representation, relative to that in vivo, of a specialized cell type and use this comparison to generate a functionally improved in vitro PC population. We predict that the generation of rationally improved cellular models will facilitate mechanistic exploration of specific disease-associated genes in their respective cell types. Keywords: Single-cell RNA-seq, Chemical biology, Stem cell-derived models, Paneth cell, Intestinal organoid, Intestinal stem cell, Differentiation, Systems biology * Correspondence: bemead@mit.edu; jmkarp@bwh.harvard.edu Benjamin E. Mead and Jose Ordovas-Montanes contributed equally to the work as co-first authors. Alexandra P. Braun and Lauren E. Levy contributed equally to the work as co-second authors. Division of Engineering in Medicine, Department of Medicine, Brigham & Women’s Hospital, Harvard Medical School, Boston, MA, USA Full list of author information is available at the end of the article © Karp et al. 2018 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Mead et al. BMC Biology (2018) 16:62 Page 2 of 24 Background and have been identified in animal models to cause clear Intestinal organoids, derived from intestinal stem cells PC phenotypes, including lower DEFA expression [20], (ISCs) and composed of ISCs, Paneth cells (PCs), enter- defects in autophagy, granule formation, and secretion oendocrine cells (EECs), goblet cells, and absorptive [19, 21], and uncompensated ER stress [22]. While in vivo enterocytes, have been invaluable to the study of intes- models currently provide the most physiologically repre- tinal biology [1]. Recent advances in massively parallel sentative system to probe PC biology, they are inherently single-cell RNA-sequencing (scRNA-seq) [2] have en- complex and poorly scaled, hindering basic research into abled the cataloging of cell types and states of the molecular mechanisms of disease and the potential for murine small intestinal epithelium [3] and intestinal scalable therapeutic screening. Recently, conventional organoids [4], offering extensive insight into tissue het- intestinal organoids were used to describe the dynamics of erogeneity, specifically within subsets of rare secretory PC degranulation in response to multiple agonists [23] cell populations. However, there have been no formal and to assess PC suppression of enteric pathogens [24]. comparisons of how the in vitro intestinal organoid While these organoid studies are arguably more repre- condition recapitulates the defined in vivo cell types. sentative than other in vitro systems, the question of While the generation of comprehensive cellular atlases physiological fidelity of this heterogeneous system re- has become a major focus of a global effort to map mains unanswered, especially given that the timescales tissues in humans, model organisms, and derived to derive conventional organoids is typically less than a organoids at single-cell resolution [5], the challenge of week, while in vivo the lifespan of a PC is on the order how to functionally investigate key insights from cell of several weeks. types in vivo, or even more simply to confirm the To improve the representation of specific cell types in high-fidelity representation of these states in existing intestinal organoids, investigators have utilized cellular model systems remains [6, 7]. engineering approaches starting with ISCs to derive Intestinal organoids are a compelling system with multiple enriched or specialized models. These include which to study specialized cells of the epithelium. They enterocytes with improved intestinal ion transport [25], are self-organizing, stem cell-derived structures, which, epithelial monolayers capable of secretion and IgA to a reasonable degree, resemble their in vivo counter- transcytosis [26], and organoids enriched for the rare part, can be rapidly grown, and are amenable to many secretory EEC population [27]. However, in each in- biochemical and genetic perturbations [8]. Recent work stance, there has been no global comparison of the has demonstrated the utility of organoids in assessing extent to which intestinal organoids, or further special- bulk phenotypes that are readily observed and easily ized derivatives, recapitulate defined in vivo cell types selected for, such as phenotypes of cystic fibrosis and the and states. Moving beyond the generation of in vivo study of cancer-associated mutational signatures [9, 10]. tissue maps towards mechanistic insights, particularly in However, the application of organoid models to the study disease settings, will require an understanding of how of complex disease, such as polygenic inflammatory the in vitro organoid models utilized for such studies disease, has been limited. In such instances, the subtler represent the cell types and states identified beyond phenotypes, such as those present in inflammatory bowel single marker genes. disease (IBD), may not manifest if the originating cell state Here,weprovideaglobalcomparisonbetween thein present in vivo is not accurately represented within an vivo cellstatesofthe murine smallintestinalepithelium organoid model. This challenge is particularly clear in IBD and the in vitro conventional intestinal organoid, and [11], where loci identified through genome wide associ- establish a systematic workflow for improving the physio- ation studies (GWAS) have proven difficult to efficiently logical representation of stem cell-derived cell states to examine through the use of in vivo animal models. enable the creation of high-fidelity in vitro models. Taking For instance, PC dysfunction is implicated in Crohn’s the PC as a test case, we utilize massively parallel single-cell disease, a subset of IBD typically afflicting the small bowel transcriptomics (scRNA-seq; Seq-Well) [28] to bench- [12]. Co-localized with LGR5 ISCs of the small intestinal mark the conventional organoid model against its in vivo crypts, long-lived PCs [13] support maintenance of the counterpart and identify differences in developmental ISC niche, producing the Wnt and Notch signaling ligands pathway signaling between in vitro and in vivo cell states. WNT3 and DLL4 [14], and are potent modulators of the Single-cell transcriptomic approaches were key in enab- gut microflora through secretion of multiple antimicro- ling this study as epithelial cell types are challenging to re- bials including lysozyme (LYZ), phospholipase A2 group liably and prospectively isolate by fluorescence-activated 1B (PLA2G1B), angiogenin ribonuclease A family member cell sorting (FACS) due to the absence of robust surface 5 (ANG5), and alpha-defensins (DEFAs), amongst others markers and the spectrum of differentiation states present. [15]. Multiple allelic variants of NOD2, ATG16L1,and This profiling guides the rational augmentation of signal- XBP1, are associated with ileal Crohn’s disease [16–19] ing pathway activity during stem cell differentiation with a Mead et al. BMC Biology (2018) 16:62 Page 3 of 24 small molecule chemical induction method we previously our study to relate organoid-derived cell states to in vivo validated to enhance global Lyz gene expression [29]. We PCs. They are fully inclusive of the 14 high confidence validate our approach by generating an enhanced in vitro markers described for Paneth cells from the terminal physiological mimic of the in vivo PC and provide a ileum in the recently published mouse small intestinal detailed characterization of the derived cell state through atlas [3]. Of note, we extended our gene list beyond truly morphologic, proteomic, transcriptomic, and functional specific marker genes that are not expressed in other cell assays based on known signatures of in vivo PCs. Further- types as we were interested in a more comprehensive set more, we use our enhanced model and findings from its of PC-enriched genes for further comparison. transcriptomic and proteomic characterization to identify We next performed scRNA-seq using Seq-Well on Nupr1 as a potential stress-response factor that facilitates conventional organoids derived from a single donor the survival of PCs, demonstrating the improved ability ISC-enriched state (Fig. 1a). Beginning with murine small to examine gene function in vitro within a more repre- intestinal crypts, we directly enriched for LGR5+ ISCs sentative cell type. over 6 days following isolation within a Matrigel scaffold and medium containing recombinant growth factors EGF Results (E), Noggin (N), and R-spondin 1 (R), small molecules Using the PC to benchmark cell type representation of CHIR99021 (C), and valproic acid (V), as well as Y-27632 conventional organoids against their in vivo counterparts for the first 2 days to inhibit rho kinase and mitigate Conventional intestinal organoids produced from the anoikis, as previously described (ENR+CV) [29]. To spontaneous differentiation of ISCs have been used to ensure reproducibility within our system and limit the risk study PCs in vitro in multiple contexts [23, 24]. These in of interference in our chemical induction approach, we vitro PCs exist as part of a heterogeneous system, yet to conducted our study exclusively with recombinant growth be rigorously benchmarked against their in vivo counter- factors and not cell line-derived conditioned media. Cells parts. To better understand the composition of PCs were passaged into conventional ENR culture for an within conventional organoids and how well those PCs additional 6 days to allow multi-lineage differentiation approximate their in vivo counterparts, we sought to and produce stem cell-derived in vitro PCs. Following globally compare the conventional organoid-derived PCs scRNA-seq, we computationally identified six clusters and their in vivo counterparts through a single-cell tran- (amongst 2513 cells × 16,198 genes meeting quality scriptomic approach (Fig. 1a). standards, see Methods) in ENR organoids, which we To relate the organoid-derived PC state to in vivo PCs, label as ENR1-4, and EEC-1 and -2 for two EEC we first generated an unbiased reference in vivo scRNA-seq types (Fig. 1d). We identified ENR-4 as the cluster data set. We performed massively parallel scRNA-seq using most enriched for Lyz1 and our PC reference gene set the recently developed Seq-Well platform [28] on epithelial (effect size 0.721, ENR-4 vs. all ENR, *t test p <2.2 × −16 cells from the ileal region of the small intestine acquired as 10 ; for effect size details see Methods)(Fig. 1e, f). two biological replicates (see Methods). We assessed quality Having identified ENR-4 as the cell state of interest metrics for the number of genes, unique molecular identi- in organoids, we directly compared the top 200 most fiers (UMIs), mitochondrial genes, and ribosomal genes, all PC-like cells in ENR-4 to in vivo PCs by performing of which fell within expectations (all cells average: 1043 differential expression analysis (Fig. 1g). In comparing genes, 2168 UMIs, 5.4% ribosomal genes, 10.4% mito- the two cell types, it became evident that the majority of chondrial genes). UMI-collapsed cell-by-gene (7667 genes enriched by in vivo PCs were defensins and antimi- cells × 17,505 genes) expression matrices were analyzed crobials, including Defa22, Defa21, Zg16, Ang4, Defa3, −37 using Seurat (see Methods), performing dimensionality and Lyz1 (all p <2.92×10 , bimodal test, Bonferroni reduction, graph-based clustering, and deriving lists of corrected for multiple comparisons) (Fig. 1g, h). ENR-4 cluster-specific genes in order to identify PCs. Within cells were enriched for Chgb, an enteroendocrine the spectrum of cell types, we identified two clusters (2 marker, and translational biosynthetic genes likely and 11) enriched for Lyz1 expression (Fig. 1b, c), of indicative of the high rates of proliferation present which we determined cluster 11 to be fully mature PCs in ENR organoids (Fig. 1g). We further note the (n = 189 cells) based on uniform expression of a set of difference in genes arising from non-sex matched associated antimicrobial peptide marker genes such as comparison, like Xist, as a limitation of our compari- Defa22, Defa21,and Ang4 (receiver operating charac- son between a single donor for organoid derivation. teristic (ROC) test, area under the curve (AUC) > 0.99 Beyond these selected genes, we note a global reduc- for markers listed; cluster 11 average: 866 genes, 3357 tion in the fraction of the transcriptome of ENR-4 UMI, 3.5% ribosomal genes, 4.8% mitochondrial genes) cells producing the total cadre of in vivo PC marker (Additional file 1: Table S1). We further utilized these genes (effect size 1.25, InVivo vs. ENR, *t test p < −16 genes (genes with AUC > 0.65 for in vivo PC) throughout 2.2 × 10 ), suggesting that the current in vitro Mead et al. BMC Biology (2018) 16:62 Page 4 of 24 ab c d ef gh Fig. 1 (See legend on next page.) Mead et al. BMC Biology (2018) 16:62 Page 5 of 24 (See figure on previous page.) Fig. 1 Transcriptional benchmarking of in vitro Paneth cells (PCs) to in vivo. a Schematic of intestinal epithelial cell isolation from terminal ileum for unbiased identification of in vivo PC signature genes, and system for intestinal stem cell (ISC) enrichment to characterize in vitro PCs, via high-throughput scRNA-seq. b Marker gene overlay for binned count-based expression level (log(scaled UMI + 1)) of Lyz1, a canonical PC marker gene, on a tSNE (t-stochastic neighbor embedding) plot of 7667 small intestinal epithelial cells isolated from the terminal ileum; receiver operating characteristic (ROC)-test area under the curve (AUC) = 0.995, n = 2 mice, independent experiments (Additional file 1:Table S1). c Violin plot for the count-based expression level (log(scaled UMI + 1)) of Lyz1 across clusters identified through shared nearest neighbor (SNN) analysis (see Methods) over small intestinal epithelial cells; n = 196 cells in cluster 11, 7667 cells in total. d A tSNE plot of 2513 cells, with clusters identified through SNN (Additional file 1: Table S1 for full gene lists with ROC-test AUC > 0.60) from conventional ENR organoids; n =6 wells of ENR organoids. e Marker gene overlay for binned count-based expression level (log(scaled UMI + 1)) of Lyz1 on a tSNE plot from; ROC-test AUC = 0.856. f Violin plot of expression contribution to a cell’s transcriptome of PC genes across ENR organoid clusters from (d) (In vivo PC gene list −16 AUC > 0.65, Additional file 1: Table S1); effect size 0.721, ENR-4 vs. all ENR, *t test p <2.2 ×10 . g Row-normalized heatmap of top differentially expressed genes using bimodal test over single-cells from the top 200 PC-like cells from ENR-4 and the 196 in vivo PCs (cluster 11, from (c)); *bimodal −16 test, all displayed genes p <1.89 ×10 or less with Bonferroni correction. h Violin plots for the count-based expression level (log(scaled UMI + 1)) of −37 Lyz1, Ang4,and Defa3 in ENR and in vivo PCs; *bimodal test, all p <2.92 ×10 or less with Bonferroni correction. i Violin plot of expression contribution −16 to a cell’s transcriptome of PC genes (effect size 1.25, InVivo vs. ENR, *t test p <2.2 × 10 ), Wnt pathway (effect size 0.559, InVivo vs. ENR, *t test −8 −7 p <2.035 ×10 ) and Notch pathway (effect size −0.500, InVivo vs. ENR, *t test p <5.25 ×10 ) genes (see Additional file 2: Table S2 for gene lists) organoid-derived PCs are suboptimal for physio- plateaued by day 6 in ENR+CD versus ENR popula- logical studies (Fig. 1i). tions. Lgr5 expression in ENR+CD at 2 days versus Modulating key developmental pathways of stem ENR showed an insignificant plateau of expression, cell-derived systems has emerged as a paradigm in which trended down by 6 days. This may be indicative bioengineering to rationally generate cell types for basic of an expansion in ‘label-retaining’ secretory precur- research and therapeutic aims [30, 31]. Specifically, modu- sors [35]. The precursor population ENR + CV had no lating Wnt and Notch signaling has been suggested in the significant difference in PC or ISC markers relative to literature to increase the frequency and magnitude of Lyz1 ENR. ThesignificantincreaseinPCgeneexpression expression and protein in ISC-derived cells [29, 32–34]. in ENR + CD relative to ENR and ENR+CV over the Leveraging the single-cell transcriptomes of our in 6-day treatment suggests rapid enrichment following vitro- and in vivo-derived PCs, we confirmed that CI, supporting our hypothesis that alterations in Wnt Wnt target genes are enriched in vivo relative to in and Notch result in superior PC enrichment in vitro. vitro PCs (effect size 0.559, InVivo vs. ENR, *t test To phenotypically describe PC enrichment following −8 p <2.035 ×10 ) and Notch target genes were decreased CI, we performed imaging and immunocytochemistry −7 (effect size −0.500, InVivo vs. ENR, *t test p <5.25×10 ) for PC-associated features. After 6 days of ENR + CD, (Fig. 1i, Additional file 2: Table S2). As a result, we sought cell populations exhibited darkened annular morphology to comprehensively test if driving Wnt and inhibiting consistent with increased numbers of granule-rich cells Notch truly results in a more physiologically representa- (Additional file 3: Figure S1A). Confocal microscopy of tive PC versus the organoid-derived PC, beyond bulk mea- whole cell clusters stained for anti-DEFA and anti-LYZ sures of increased Lyz1 expression. showed an increase in LYZ+ and DEFA+ cells in ENR + CD compared to both ENR and ENR + CV (Fig. 2c). Rationally guided chemical induction of Wnt and Single-cell counting of confocal imaging showed a inhibition of Notch drives PC marker enrichment significant increase of DEFA and LYZ co-staining cells Beginning with an LGR5+ ISC-enriched population (ENR in ENR + CD (20–30% of cells) versus either ENR or +CV), we sought to profile how the modulation of Wnt ENR + CV (both < 5%; adj. p = 0.0001) (Additional file 3: and Notch signaling through small molecule inhibitors Figure S1B). Additionally, normalized z-axis profiles of in- would alter the in vitro PC state, as suggested by our dividual co-staining cells within cell clusters revealed a con- transcriptional profiling. We performed chemical induction sistent distribution of DEFA (luminally polarized) and LYZ (CI) using the previously identified compounds C to drive (diffuse) (Additional file 3:FigureS1C 1–3). High-resolution Wnt signaling and DAPT (D), a gamma-secretase inhibitor, fluorescent imaging of individual co-staining cells from to inhibit Notch (ENR+CD) (Fig. 2a)and measured bulk freshly isolated small intestinal crypts (in vivo equivalent) gene expression of PC (Lyz1, Defa1, Mmp7)and ISC(Lgr5) and 6-day ENR + CD-treated cells showed a similar markers every 2 days for a total of 6 days (Fig. 2b). ENR polarized distribution of LYZ- and DEFA-stained gran- +CD-treated cells had statistically significant increases ules, although freshly isolated cells appeared to be in Lyz1 (adj. p = 0.005, see Methods)and Mmp7 (adj. more granular than CI-PCs (Fig. 2d). p = 0.005) within 2 days compared to ENR, with differ- To confirm the extent of enrichment seen in whole ences plateauing around 4 days. Defa1 (adj. p = 0.004) population imaging, the prevalence of PCs in ENR + CD expression was significantly increased by day 4 and relative to ENR was assessed by flow cytometry over the Mead et al. BMC Biology (2018) 16:62 Page 6 of 24 ab c de f gh Fig. 2 (See legend on next page.) Mead et al. BMC Biology (2018) 16:62 Page 7 of 24 (See figure on previous page.) Fig. 2 Establishing chemically induced Paneth cell (PC)-enriched cultures. a Schematic of small molecule-driven differentiation of LGR5+ ISCs (C - CHIR99021, D - DAPT) and non-specific differentiation. b mRNA expression of PC (Lyz1, Defa1, Mmp7) and ISC (Lgr5) markers relative to ENR, for ENR+CV and ENR + CD at 2 (D2), 4 (D4), and 6 days (D6) (n = 3 biological replicates; two-way ANOVA with multiple comparison test vs. ENR; ** adj. p < 0.01, *** adj. p < 0.001). c Representative confocal imaging of whole cell clusters for PC antimicrobials following 6 days in ENR+CD versus ENR and ENR+CV: stained for anti-DEFA, anti-LYZ and counterstained with DAPI and for actin (phalloidin). d High-resolution fluorescent imaging of in vivo and in vitro single cells from 6-day culture in ENR + CD shows similar morphology and antimicrobial expression: stained for DEFA and LYZ, and counterstained with DAPI and for actin (phalloidin). e Viable cell populations from ENR, ENR+CD, and ENR+CV precursor culture have distinct populations based on CD24 and LYZ content, indicative of PC maturity (n = 3 biological replicates; ENR + CV, days 4, 6, 12, n = 2 biological replicates day 8). f Volcano plot of differentially regulated proteins between 6-day (6D) ENR + CD and ENR cells shows clear enrichment in secreted and PC-associated proteins (labeled). Cut-offs are 2 standard deviations outside the mean expression level of the set and FDR < 0.05. g Rank-order log fold change of detected PC antimicrobial proteins and between 6-day ENR + CD and ENR cultures (n =4). h Rank-order log fold change of detected secretory proteins associated with EEC and goblet lineages in ENR + CD relative to ENR cultures (n =4) course of 12 days, a longer term culture than typical for samples and analyzed them in a single 10-plex by conventional organoids. We identified an in vivo PC LC-MS/MS (Additional file 4: Figure S2A). We identi- phenotype as CD24 and LYZ co-positive cells, as per fied 8015 unique proteins within all samples; each previous reports [36], and noted the presence of replicate pair (ENR + CD/ENR) was normally distributed single-positive LYZ+ or single-positive CD24+ popula- (not shown) and correlated with all others, indicating con- tions, indicative of alternative cell differentiation, imma- sistentproteomeenrichment(Additional file 4:FigureS2B). ture, or non-physiological PCs (representative populations We looked at the sample pairs in aggregate and classified Additional file 3:FigureS1D,representative gating proteins significantly enriched in ENR + CD and ENR by a Additional file 3: Figure S1E). ENR + CD had substantial false discovery rate (FDR) < 0.05 and log fold change (± 2σ) enrichment at all time points for double-positive, and (Fig. 2f and Additional file 5: Table S3). There were 249 single-positive LYZ+ or CD24+ populations relative to ENR + CD-enriched proteins, 212 ENR-enriched proteins, ENR, as well as a consistent decrease in the double and 7553 shared proteins. Known PC markers, including negative population in agreement with the PC phenotype LYZ, DEFAs, and other secretory pathway components, (Fig. 2e). Notably, both ENR and ENR + CD experience were identified as significantly enriched in ENR + CD ver- declines in total cell viability, with ENR + CD having sus ENR alone. Of known antimicrobial proteins produced greater survival at longer times, suggesting both a reduc- by PCs, we detected 10 DEFAs, 5 CRS peptides, 6 ribonu- tion in anoikis, a potentially physiological ‘long-lived’ PC cleases, 12 lectins, LYZ1, and PLA2G1B with differential phenotype in ENR + CD versus ENR, or an enhance- abundance between ENR + CD and ENR (Fig. 2g). Each ment in niche-supporting functionality (Additional class of antimicrobials had at least one ENR + CD enriched file 3: Figure S1F). Overall, imaging and flow cytometry protein (+ 2σ), with the ribonucleases significantly enriched demonstrate a significant increase in cells morphologically and a majority of the lectins and DEFAs unregulated be- resembling in vivo PCs with respect to granularity, polar- tween the two conditions. Proteins associated with the EEC ity, and antimicrobial co-expression in ENR + CD com- lineage (secretogranins, chromogranins, and neuropeptides) pared to conventional ENR organoids (Fig. 2c–e and were also enriched in ENR + CD, in addition to multiple Additional file 3:FigureS1A–F). other secreted components, including Wnt ligands, and the complement pathway components C3 and CFI (Fig. 2h). In Chemically induced PC proteome is enriched for sum, we see a broad diversity of PC-associated antimicro- components of secretory lineages bials with some enrichment of EEC-associated proteins in With ENR + CD apparently providing a more prevalent ENR + CD relative to ENR. and physiological PC population, we sought to more Additionally, we characterized enriched biological func- globally characterize the differences between in vitro tions, cellular compartments, and molecular functions PCs (ENR vs. ENR + CD) at 6 days. Because our single using DAVID v6.8 and the gene ontology database. All cell transcriptomic comparison revealed that many of sets had high database coverage (greater than 85%) of the differential genes between PCs in conventional orga- queried proteins. The ENR + CD proteome is signifi- noids and in vivo were lineage-defining protein prod- cantly enriched for extracellular and protein process- ucts, we sought to assess the total intracellular proteome ing compartments and secretory-associated functions between the conventional organoid and our chemically (Additional file 4: Figure S2C), while the ENR prote- induced model through liquid chromatography mass ome favors translation, intracellular compartments, spectrometry (LC-MS/MS)-based proteomics. We quan- and translational activities (Additional file 4:Figure tified relative protein abundance using isobaric mass S2D). Of note, there are the extracellular exosome and tag labeling from four ENR and four ENR + CD calcium ion-binding associated proteins in the ENR + Mead et al. BMC Biology (2018) 16:62 Page 8 of 24 CD proteome that are indicative of the intestinal cells expressing antimicrobial Lyz1, Defa24, Defa3, Mmp7, epithelial secretory phenotype (for a complete list of and EEC marker Chga, and that ENR enriched for absorp- DAVID enrichments, refer to Additional file 6:Table tive marker Fabp2-expressing cells (Fig. 3b). S4). These functional enrichments further support the To assess subpopulation structure and provide a notion that the ENR + CD-cultured organoids are more robust measure of composition beyond canon- enriched in secretory cells, including PCs, although it ical marker genes, we performed unsupervised KNN does not rule out potential co-enrichment for the EEC graph-based clustering on the captured cells (Fig. 3c, d lineage. Finally, we sought to identify transcription and Additional file 1: Table S1 for full gene lists), factors (TFs) that may mediate PC-specific differenti- distinguishing four clusters in each treatment condi- ation using GSEA [37, 38]withthe MSigDB transcrip- tion. We then scored individual clusters according to tion factor target (v5.2) gene set database [39]witha the amount of the transcriptome within each cell dedi- moderately conservative cutoff (see Methods). We cated to synthesizing the respective enriched proteins generated an enrichment map [40, 41]ofseveral TF from the bulk proteome data. We observed that ENR targets significantly enriched in both the ENR + CD + CD clusters yield a significant enrichment for those and ENR proteomes. In ENR + CD, the nuclear recep- proteins detected in the up-regulated proteome (effect −16 tors for progesterone, aldosterone, and glucocorticoid, size 1.38 ENR + CD vs. ENR clusters, p <2.2 ×10 ) as well as the cellular differentiation-implicated TAL1, and that the down-regulated proteins were enriched in RP58, and NRSF, were significantly enriched. In ENR, the ENR and ENR + CV conditions (Fig. 3d, e and data the primary known enrichment was for the cell cycle and not shown). Intriguingly, at the level of clusters, the proliferation-related E2F TF family (Additional file 4: upregulated proteome was not evenly distributed Figure S2E). These potential TFs are consistent with across all cells in ENR + CD, but was most enriched in CI-PC treatment driving expected terminal differentiation cluster ENR + CD-4, which represented approximately of specialized cells, as opposed to conventional organoid 10% of ENR + CD cells (effect size 2.40 ENR + CD-4 −16 culture, which supports a broad mix of intestinal epithelial vs. all cells, p <2.2 ×10 )(Fig. 3d, e). cells, including proliferating populations. Furthermore, To address ENR + CD composition and how it relates to this analysis suggests potential targets, such as progester- conventional organoids, we interrogated the expression of one, aldosterone, and glucocorticoid, to modulate the Lyz1, Chga, and other selected genes across each cluster differentiation programs of this secretory cell population (Fig. 4a). We noted that clusters ENR-4 and ENR + CD-4 in future studies. shared expression of Lyz1, Defa24, Defa3,and Mmp7, yet ENR + CD-4 cells produced significantly more of each −74 Single-cell RNA sequencing profiles heterogeneity of canonical PC gene (bimodal test, p <6.80×10 for genes chemically induced PCs, revealing subsets with improved listed, Bonferroni corrected for multiple comparisons). transcriptional similarity to in vivo Furthermore, both ENR-4 and ENR + CD-4 cells lacked With the apparent co-enrichment of canonical PC and expression of EEC genes like Chga, which was observed in EEC proteins in the ENR + CD proteome, we sought to the EEC-1 and EEC-2 clusters arising from mixed-grouping identify whetherweproduceahomogenouspopulationof of the sample, as well as in ENR + CD-2 and ENR + CD-3 mixed-lineage secretory cells or a spectrum of unique cell (Fig. 4a). Altogether, this suggests that ENR + CD drives states between EEC and PC. We performed scRNA-seq PC differentiation while also inducing a secretory tran- using the Seq-Well platform on cells from ENR + CD and sition state (ENR + CD-2 and 3) expressing a mix of PC the precursor ENR + CV conditions to analyze alongside and EEC marker genes (Additional file 1:Table S1 for conventional ENR organoids. To ensure experimental full gene lists). robustness, we assessed quality metrics for the number of We next sought to compare the states generated in genes, UMIs, mitochondrial genes, and ribosomal genes by vitro to those observed in vivo with our refined system. cluster, all of which fell within expectation (Additional file 7: Using the gene list of in vivo PC markers and further de- Figure S3). UMI-collapsed digital gene expression matrices fining a list for in vivo EECs (see Methods) captured on were analyzed using Seurat (see Methods), and displaying the Seq-Well platform (Additional file 1: Table S1), we all three treatments (ENR + CV, ENR, ENR + CD) in tSNE observed that the percentage of a cell’stranscriptome space demonstrated clear separation between each condi- dedicated to synthesizing defining Paneth genes was tion (Fig. 3a), illustrating that the unique transcriptional significantly enriched relative to ENR-4 in clusters ENR −5 differences induced by each treatment were conserved + CD-2, -3, and -4 (effect size 0.15, p <3.43 × 10 ;ef- −16 across all cells. Plotting key genes demonstrated that, as fect size 0.829, p <2.2 ×10 ; effect size 2.52, p <2.2 × −16 expected, all cells expressed high levels of Epcam,that 10 , respectively) with an increase in expression of ENR + CV cells had enhanced Mki67, a marker of prolifer- EEC genes across ENR + CD-1, -2, and -3 but not ENR –16 ation, that the ENR + CD condition led to enrichment of + CD-4 (effect size 1.30, p <2.2 ×10 ; effect size 1.82, Mead et al. BMC Biology (2018) 16:62 Page 9 of 24 a b cd e Fig. 3 Single-cell RNA-sequencing reveals cellular composition across treatments and origins of proteomic data. a A tSNE plot of single cells derived from ENR + CV (n =985 cells), ENR (n =2544 cells), and ENR +CD (n = 2382 cells) harvested at day 6 of differentiation, colored by treatment; n =6 wells for each condition. b Marker gene overlays (on plot from (a)) for binned count-based expression level (log(scaled UMI + 1)) of individual genes of interest. c A tSNE plot, with clusters identified through SNN graph-based clustering (see Additional file 1: Table S1 for marker gene lists), highlighting distinct cell states within each organoid; opacity of density clouds correspond to the Paneth cell score of ENR-4, ENR + CD-3, and ENR + CD-4 clusters (see Fig. 4b). d Violin plot of expression contribution to a cell’s transcriptome of ENR + CD proteome-enriched genes across organoid clusters from (c) (see Additional −16 file 1:Table S1 forfullgene list);effect size 2.40 ENR+CD-4 vs.all cells, p <2.2 ×10 . e Frequency of each cluster observed within each organoid condition as a fraction of the total cells in each condition –16 −16 p <2.2 ×10 ; effect size 1.118, p <2.2 ×10 ; effect size ENR + CD-4 cells relative to in vivo PCs demonstrated a 0.0465, p = 0.2339, respectively) (Fig. 4b). Notably, ENR + striking similarity relative to the difference observed be- CD-4 cells (~10%) had a three-fold increase in the tran- tween in vivo and ENR-4 cells (PC fraction of in vivo scriptional resemblance to in vivo PCs relative to ENR-4 transcriptome: effect size 0.237 InVivo vs. ENR + CD-4, (53.4% of transcriptome ENR + CD-4 vs. 16.5% of tran- p < 0.0055; effect size 1.25 InVivo vs. ENR-4, p <2.2 × −16 scriptome ENR-4) (quantification of Fig. 4b). Furthermore, 10 (Additional file 1:Table S1). 45% of ENR + CD cells express a secretory PC-like tran- In Fig. 4c, we present a heatmap of scaled expression scriptional phenotype that is at least two-fold enhanced values for the top genes (AUC > 0.65) used for the in relative to conventional organoids (33.9% of transcriptome vivo Paneth score across ENR-4, ENR + CD-4, and the in ENR + CD-3 and -4 vs. 16.5% ENR-4). Comparing the vivo cluster used to define PCs. We observed that the Mead et al. BMC Biology (2018) 16:62 Page 10 of 24 Fig. 4 (See legend on next page.) Mead et al. BMC Biology (2018) 16:62 Page 11 of 24 (See figure on previous page.) Fig. 4 Transcriptional identity of chemically induced Paneth cells (CI-PCs) within conditions and related to in vivo PCs. a Violin plots for the count-based −74 expression level (log(scaled UMI + 1)) of selected genes across called clusters, colors correspond to clusters in Fig. 4c;*t test, p <6.80× 10 or less with Bonferroni correction, for Lyz1, Defa24, Defa3,and Mmp7 ENR + CD-4 relative to ENR-4. b Violin plot of expression contribution to a cell’s transcriptome of in vivo PC and enteroendocrine marker-cell genes (see Additional file 1: Table S1 for full gene list, AUC > 0.65); effect size 2.52 ENR + CD-4 vs. −16 ENR-4, p <2.2 × 10 for PC score; effect size 0.0465, p = 0.2339 ENR + CD-4 vs. ENR-4 for enteroendocrine cell score. c Row-clustered heatmap of z-scores (−2.5 to 2.5; purple to yellow) for defining genes (n = 69 with AUC > 0.65 of in vivo PCs, see Additional file 1:Table S1 forfullgene list) across top 200 cells for PC score (Fig. 5b) from ENR-4 and ENR + CD-4 conditions compared to two biological replicates of in vivo PCs from the terminal ileum (n = 196 cells) enhanced PC phenotype in ENR + CD-4 (effect size around 6 hours post-wash (two-way ANOVA, stimulant −16 1.144 ENR + CD-4 vs. ENR-4, p < 2.2 × 10 ) correlated p < 0.0001, time-point p < 0.0001) (Fig. 5b). The ob- with a greater expression of signature genes, such as served PC secretion in response to CCh is consistent Lyz1, Lyz2,and Defa5, and a greater diversity of anti- with observations made in ex vivo crypts, though over microbial peptide genes, such as Ang4, Defa3, and the appreciably longer time scales, likely due to the added metalloprotease Mmp7. diffusion barrier of the organoid structure and matrigel To confirm and extend our findings of pathway-based [44]. We next identified how LYZ secretion changes modulation, we scored clusters for enrichment or depletion over the course of differentiation. Beginning with an of canonical growth factor-induced pathways. CHIR99021 ISC-enriched population, we assayed for secreted LYZ activates the Wnt pathway, and we observed a significant in cell culture supernatants every 2 days for 6 days of enrichment for Wnt target genes in all CI-PC clusters ENR + CD culture, following a 24-h stimulation with −16 (effect size > 0.999, p < 2.2 × 10 for all ENR + CD clus- CCh or without (basal collection/non-stimulated). ters vs. ENR-4) (Additional file 8: Figure S4A). While Notable increases in functional secretion (stimulated DAPT is a Notch pathway inhibitor, levels of Notch target relative to basal) occurred at days 4 and 6 (two-way genes were largely greater than or equivalent to ENR-4 ANOVA, stimulant p < 0.0001, time-point p < 0.0001) cells across CI-PC clusters, except for significant depletion (Fig. 5a). Compared to conventional organoids and −16 in ENR + CD-4 (effect size −0.658, p <2.2 ×10 ENR + ISC-enriched precursors, ENR + CD secreted significantly CD-4 vs. ENR-4) (Additional file 8: Figure S4B). This more basal LYZ (p < 0.0001) and was the only population suggests that complete Notch suppression is key for PC that showed grossly measurable CCh-induced secretion differentiation distinct from an EEC fate. Additionally, (adj. p = 0.03) (Fig. 5c). This result is consistent with the given the recognized role for distinct respiratory potential observed enrichment between chemically induced popula- in enterocytes, ISCs, and PCs, we scored cells across tions relative to conventional. respiratory electron transport genes [42, 43]. ENR + CD-4 Based on the broad spectrum of antimicrobials de- had the lowest cluster score relative to all cell subsets tected proteomically, transcriptionally, and functionally, –16 (effect size −1.4649, p <2.2 ×10 )(Additional file 8: we hypothesized that ENR + CD possess greater bacteri- Figure S4C). Together, this suggests that Wnt signaling cidal effects than conventional organoids. We assayed is necessary but not sufficient to specify the mature PC for bacterial growth modulation by suspending cell clus- phenotype and that Notch and metabolic conditions ters with common laboratory strains of gram-negative play a larger role in the decision between PC and EEC and gram-positive bacteria in exponential growth. fates. CI-PCs significantly suppressed growth of gram-positive L. lactis MG1363 (adj. p = 0.0001), which did not occur Chemically induced PCs mimic in vivo stimulant-induced with conventional organoids, indicative of increased secretion and demonstrate selective modulation of PC-associated antimicrobial activity (Fig. 5d). Both ENR bacteria in co-culture (adj. p = 0.0005) and ENR + CD (adj. p = 0.01) co-culture In addition to our morphological, proteomic, and tran- showed significant increase in gram-negative E. coli scriptional characterization of PC phenotype in ENR + MG1655 growth but no appreciable effect was observed CD and ENR, we sought to measure physiological on the growth of gram-positive E. faecalis V583 versus function by assessing stimulant-induced secretion of bacteria alone (Fig. 5d). While this assay simplifies the antimicrobials. We assessed the dynamics of LYZ PCs’ physiological environment and may not be a direct accumulation in the supernatant media of cultures proxy for strain-specific growth modulation, it does following media wash, basally and after stimulation demonstrate that the PC-enrichment of ENR + CD ver- with carbachol (CCh), a cholinergic agonist known to sus conventional organoids enables detectable in vitro induce PC secretion [44]. CCh (10 μM) induced a rapid bacteria species-specific PC antimicrobial response, accumulation of LYZ within 2 hours that plateaued opening avenues for future experimentation. Mead et al. BMC Biology (2018) 16:62 Page 12 of 24 a bc de Fig. 5 Chemically induced Paneth cells (CI-PCs) are functional in response to host and microbial stimuli. a Supernatant LYZ from 24-h basal and 10 μm CCh-stimulated LYZ cells at varying number of days in ENR + CD culture (top). DNA content from matched samples (bottom) (n = 8 well replicates; SEM error bars too small to visualize). b Supernatant LYZ from 6-day ENR + CD collected basally and following 10 μm CCh stimulation for 0.5, 2, 4, 6, and 24 h (top). DNA content from matched samples basally and following 10 μm CCh stimulation (bottom) (n = 8 well replicates). c 24-h basal (non-stimulated) and 10 μm CCh-stimulated LYZ secretion in 6-day ENR + CD versus ENR and ENR + CV (n = 8 well replicates; two- way ANOVA with multiple comparison test; ns non-significant, * adj. p < 0.05, **** adj. p < 0.0001). d 4-h co-culture of freshly passaged 6-day ENR and ENR + CD cells and select gram-negative and gram-positive aerobic bacteria (n = 13 well replicates; two-way ANOVA with multiple comparison test, * adj. p < 0.05, *** adj. p < 0.001, **** adj. p < 0.0001). e Normalized cellular viability, caspase activity per viable cell, and cytotoxicity per viable cell from 24-h and 48-h ENR and ENR + CD co-cultures at specified mixing ratios (n = 9 well replicates from three biological donors; one sample t test,* p <0.05, ** p <0.01, *** p <0.001, **** p <0.0001) Chemically induced PCs provide niche support and co-culture viability, caspase activity, and cytotoxicity 24 enhance conventional organoid survival and 48 h following re-plating in ENR media. If there was Beyond the generation of antimicrobial peptides, PCs no appreciable interaction, positive or negative, between provide niche support for ISCs. We sought to test if the two populations we would expect to see a linear CI-PCs provided niche factors known to drive epithelial trend of every measured variable throughout mixing regenerative turnover. We performed co-culture experi- ratios. However, we observe a significant positive inter- ments, mixing and re-plating cell populations derived action where the presence of both populations drives an from 6 days of ENR or ENR + CD culture and assayed overall increase in cellular viability, beginning at 24 h Mead et al. BMC Biology (2018) 16:62 Page 13 of 24 (one-sample t test 1:1 p = 0.037) and increasing at 48 h for 2 days following a 6-day course of ENR + CD, where (one-sample t test 1:1 p =0.001 and 1:3 p < 0.001) again a profound, but not total, decline in cellular viability (Fig. 5e). This is likely due to a significant decrease in was observed. Further, it appears that TFP treatment is se- overall apoptosis relative to the total cell population lectively more toxic to PC and PC-progenitor populations (one-sample t test 24 h 1:1 p = 0.004 and 1:3 p = 0.032, relative to non-PC populations (Fig. 6d). In total, this ini- 48 h 1:3 p = 0.003), and unrelated to changes in cellular tial investigation suggests that NUPR1 may be a critical cytotoxicity. We believe that the presence of a TF in PC development and survival, which carries PC-enriched population (from ENR + CD) is driving therapeutic implications which we will seek to validate this effect by providing increased soluble regenerative with expanded gene-perturbation studies in vitro and in factors to the ISC population in ENR organoids, in- vivo in future work. creasing the generation of new cells, and resulting in a lower overall rate of apoptosis. Discussion We sought to directly compare a specific cell type present Mapping of in vivo PC-associated transcription factors to in vivo to that derived in an intestinal organoid in vitro, in vitro proteome and transcriptome reveals Nupr1 as with the main goal of understanding the nature and extent important in epithelial survival of divergence between the in vitro and in vivo conditions. Finally, we sought to use this physiologically improved in Empowered by recent advances in massively-parallel vitro PC system (ENR + CD) to identify novel factors po- scRNA-seq, we employed a generalizable approach to tentially supportive of PC survival or differentiation. Using define the relation between in vivo and organoid-derived our in vivo PC and EEC gene lists, and filtering for only in vitro PCs and employ a rationally identified interven- TFs (using TFdb, downloaded September 2017) [45], we tion to improve in vitro representation through chemical identified a set of PC- or EEC-specific TFs. We mapped modulation of developmental pathways. Our scRNA-seq these TFs to our in vitro proteome (Fig. 6a and Additional approach enabled the identification of populations of file 5: Table S3), which revealed the previously unreported interest both in vivo and in vitro and the analysis of differ- NUPR1 as the most enriched PC-specific TF in ENR + ences between subpopulations which would have other- CD. This finding was supported by differential expression wise been greatly obscured in a bulk analysis of this between ENR + CD-2 (most EEC-like cells) and ENR + heterogeneous system. While others have recently used −37 CD-4 (p <3.12×10 , bimodal test, Bonferroni corrected scRNA-seq to profile the heterogeneity of intestinal orga- for multiple comparisons) (Fig. 6b). We further identified noids [4] and the murine small intestinal epithelium [3], Nupr1 in our in vivo PC populations, which showed spe- we provide a direct comparison between the model and cific and enriched expression of Nupr1 by in vivo PCs system. This comparison is key to understanding what (ROC test, AUC = 0.833) (Fig. 6b). Intriguingly, Nupr1 is a complex models, such as intestinal organoids, really repre- stress-response gene, known to promote cellular survival sent, how they may best be utilized and rational strategies and senescence through mediation of autophagy, and has for improvement. primarily been studied in the context of cancer [46–48]. In our comparison, we identified that the PC-state of Autophagy and stress response have repeatedly been im- conventional intestinal organoids poorly represents the plicated through GWAS study in PCs in IBD; however, extent and breadth of antimicrobial gene expression, and Nupr1 has only ever been reported in a single IBD GWAS that modulation of Wnt and Notch during differentiation study, and its role in PC biology has not been formally in- may improve physiological representation. Using a com- vestigated [49]. With our model, we sought to test the role bination of small molecule promotion of Wnt and inhib- of NUPR1 on in vitro PC survival through the small ition of Notch signaling that had previously been shown to molecule inhibition of NUPR1 with trifluoperazine (TFP) improve the bulk expression of PC-marker gene Lyz1 [29], [50, 51]. While genetic perturbation may provide for more we drove a secretory differentiation program and enriched specific effect measurement, we chose to use TFP as a for mature PCs with greater diversity and expression of simple, albeit less specific, modulator, as the complexity antimicrobial peptides relative to existing in vitro models involved in selecting for a genetically perturbed popula- and, thus, are more representative of in vivo PCs. Imaging tion of PCs in organoids, if Nupr1 is a survival gene, is be- of this population revealed that they are positive for the yond the scope of the present study. We first tested how antimicrobials LYZ and DEFA, clearly polarized, and different dosages impact PC differentiation in combin- granulerich, suggestive of amaturePC. This population is ation with ENR + CD for 6 days, where doses above 1 μM approximately six-fold more abundant in ENR + CD than led to near total cell death, and where the few surviving an ENR organoid, as confirmed through image quantifica- cells were primarily non-PC (Fig. 6c). This suggests that tion, flow cytometry and scRNA-seq. We further character- Nupr1 is likely critical to cellular survival during the CI ized the subpopulation enrichments of our ENR + CD differentiation process. We also tested the addition of TFP culture and directly compared it to conventional organoids. Mead et al. BMC Biology (2018) 16:62 Page 14 of 24 ab cd Fig. 6 CI-PCs reveal putative function of Nupr1 transcription factor in Paneth cell (PC) survival. a ENR + CD is enriched for in vivo PC and EEC transcription factors, including Nupr1 (n = 4). b Violin plots for the count-based expression level (log(scaled UMI + 1)) of Nupr1 across in vivo and in vitro called clusters. c Nupr1 inhibition with trifluorperazine (TFP) treatment concurrent with 6-day ENR + CD differentiation reveals dose- dependent toxicity, with preference to PCs (CD24+ and LYZ+) and PC-like (CD24+ and LYZ+) populations as assessed by flow cytometry (n =3 biological replicates). d 2-day TFP treatment following 6-day ENR + CD differentiation reveals dose-dependent toxicity, with preference to PCs (CD24+ and LYZ+) and PC-like (CD24+ and LYZ+) populations as assessed by flow cytometry (n = 5 biological replicates) We identified two subpopulations in scRNA-seq (ENR + identify as the approximately 5% of single-staining LYZ+ CD-3 and ENR + CD-4) that account for approximately cells present in ENR organoids as assessed by flow cytome- half of the ENR + CD-treated cells with a high-degree of try (Fig. 2e). This finding makes sense, given that conven- transcriptional similarity to in vivo PCs, a greater percent- tional organoid culture often occurs on the time scale of a age/matching than the ENR-subpopulation that most re- week, while in vivo PCs are relatively long lived (several sembles an in vivo PC (ENR-4). From this analysis, we weeks), develop in non-sterile conditions, and presumably believe that in vitro PCs characterized in the past [23, 24] would require longer periods of culture to reach maturity likely represent secretory precursor populations lacking the in vitro. Indeed, our ENR + CD cultures show increasing full phenotypic repertoire of the in vivo PC, which we PC populations up to 12 days of culture and may likely Mead et al. BMC Biology (2018) 16:62 Page 15 of 24 continue gaining at longer time points. While our ap- inhibition is important for mature PC development, pos- proach moves us much closer to generating the in vivo PC sibly as a balance between differentiation and cell survival. (fraction of in vivo transcriptome: effect size 0.237 InVivo Interestingly, we also see a notable gradient in cellular vs. ENR + CD-4, p < 0.0055; effect size 1.25 InVivo vs. respiration across subpopulations, lowest in the PC and −16 ENR, p <2.2 ×10 ), we still do not capture the total highest in the stem cell and EEC lineages, in agreement amount of antimicrobial peptides present in vivo, and with recent work on the metabolic differences within the propose pathways to modulate in future studies. stem cell niche [43], as another potential cue to further Evidence suggests that PC antimicrobial expression specify PC differentiation. Future studies should incorpor- and function are influenced by genetic background [52] ate temporal aspects to growth factor delivery akin to and implicated in intestinal disease, including IBD [53]. what has been shown for degradable matrices [58]to How genetic background may influence differentiation enhance purity and yield, and explore the role of meta- through this protocol is yet to be studied but especially bolic utilization in addition to growth factor signaling in prudent, as we demonstrated the ability to detect a cellular fate determination. Indeed, while we implemented broad spectrum of antimicrobial proteins and peptides an approach of chemical induction to drive model im- and their differential abundance within a PC-enriched provement, there are many other approaches which may population. Interestingly, we identified that the same be implemented to seek similar shifts in the composition subpopulation (ENR + CD-4) with the most transcrip- of organoid systems, such as using novel ligands or tuning tional overlap to our bulk ENR + CD-enriched proteome of the material supports. Further, we appreciate that many also most closely resembles the in vivo PC. While this investigators have begun using R-spondin or Wnt condi- subpopulation does not account for the majority of ENR tioned medias in their intestinal organoid cultures; how- + CD-cultured cells, it appears that ENR + CD-4 consist- ever, to ensure consistency and control of the system ently drives the PC phenotype in vitro. In addition to upon chemical induction we chose to use recombinant assessing the role of genetic background or disease state growth factors in lieu of a less-characterized media prod- on antimicrobial content, our platform also affords the uct, but this may not preclude their usage. Overall, our ability to interrogate how alterations in protein process- analyses of single cell heterogeneity show that our sys- ing and storage in PCs affects the proteome, which has tem is well positioned to further investigate the effects been shown to drive shifts in the microbiome and may of both known and unknown physiological cues on PC be implicated in disease [54, 55]. Finally, while we differentiation and function. demonstrate an enriched phenotypic spectrum of anti- One of the most important features we established with microbials and Wnt ligands, we also identified several our CI-PCs was the ability to measure PC functional neuropeptides and hormone products associated with enrichment through simple soluble assays. We demon- the EEC lineage within our system. Given that multiple strated sufficient functional enrichment in PCs such that studies have linked the differentiation of PCs and EECs enzymatic activity assays can detect stimulant-induced through a common progenitor population [56], it is secretion of antimicrobials as well as the promotion of reasonable to expect enrichment in one population the ISC niche. Moreover, microbe co-culture assays would also allow for some overlap with the other, as we with our enriched cells produce measurable and select- see in our scRNA-seq data. Future experiments will seek ive microbial growth modulation not observed using to leverage these transition states to more formally iden- conventional organoids. Co-culture strains were chosen tify genetic programs that underlie common or unique to demonstrate proof of concept of selective antimicro- developmental trajectories. bial action and assess functionality compared to con- To understand how our chemical induction led to ventional organoids. Given the results showing selective distinct secretory subpopulations within the CI-PCs, we modulation of bacterial growth, we believe that our sys- mapped Wnt, Notch, and metabolic gene sets onto each tem could serve as a tool to further probe host-microbe subpopulation. In our system, Notch-signature is highest interaction in vitro. Furthermore, it would allow for in- in the stem cells and EECs, lower in enterocytes, and vestigations of both microbial mechanisms that elicit lowest in PCs. Our system’s Wnt signature is relatively PC response (e.g., TLRs) and the properties of complex decreased in enterocytes (ENR largely) and increased in mixtures of secreted components, including multiple PCs and EECs, both of which occur predominantly in the antimicrobial proteins. Wnt-driven condition ENR + CD (CI-PCs). In total, this suggests that Wnt is necessary for ISCs to commit to PC Conclusions and EEC lineages and that future experimentation with The generation of comprehensive cellular atlases from specific synthetic Wnt ligands [57] may prove fruitful in humans and model organisms will revolutionize our distinguishing Wnt target genes that discriminatorily yield understanding of complex tissues [3]. Intestinal orga- PCs or EECs. Additionally, it is clear that strong Notch noids have already proven their value in studying human Mead et al. BMC Biology (2018) 16:62 Page 16 of 24 and murine epithelial biology. However, to rigorously Sigma-Aldrich) supplemented with growth factors EGF (E; test hypotheses of basic biological or disease mechanism, 50 ng/mL, Life Technologies), Noggin (N; 100 ng/mL, it will be essential to have simple and reliable protocols PeproTech), and R-spondin 1 (R; 500 ng/mL, PeproTech) for the generation of specialized subsets of cells which and small molecules CHIR99021 (C; 3 μM, LC Laborator- cannot be readily isolated from tissue. The representa- ies) and valproic acid (V; 1 mM, Sigma-Aldrich) was added tiveness of cell states present in organoids and the to each well. ROCK inhibitor Y-27632 (Y; 10 μM, R&D specialized cell types present in vivo [3] is an outstand- Systems) was added for the first 2 days of culture. Cells ing question with implications in mucosal immunology, were cultured at 37 °C with 5% CO , and cell culture developmental biology, and translational medicine. Our medium was changed every other day. After 6 days of single-cell genomics approach provides compelling evi- culture, crypt organoids were isolated from Matrigel by dence that organoid-derived cell populations must be mechanical dissociation. Isolated organoids were resus- validated to ensure physiological relevance, and add- pended in TrypLE Express (Life Tech) to dissociate into sin- itionally provides a rational framework for identifying gle cells, then replated in Matrigel with ENR + CV + Y cell states and their potential upstream drivers to modu- media for 2 days. Cells were once again passaged, either late cellular composition. This approach could enable ad- into freezing media (Life Tech) for cryopreservation or vances beyond conventional organoid systems to provide replated at approximately 200 organoids per well (24-well an enriched highly specialized cell population that recapit- plate) for ISC-enriched organoid expansion. ISC-enriched ulates important physiological functions of the intestinal organoids were passaged or differentiated every 6 days in epithelium and could represent an improvement in in the ENR + CV condition. To differentiate, cells were pas- vitro PC culture for the purposes of high-throughput saged as previously described, and crypt culture medium screening, the study of host-microbe interactions, bio- containing growth factors ENR only or ENR + CD (DAPT, engineering (e.g., precision gene editing), and the identifi- 10 μM; Sigma-Aldrich) was added to each well. cation of novel genetic candidates in PC function (e.g., Nupr1). With this framework, we illustrate the power and RNA extraction and qRT-PCR importance of rigorously characterizing the specialized Organoids were isolated from Matrigel in 24-well plates cell types derived in organoids to those defined in ‘atlas-le- following culture as previously described, and pellets vel’ surveys of the intestinal epithelium. were lysed in TRI reagent with RNA extracted according to the manufacturer’s protocol (T9424, Sigma). Resulting Methods RNA pellets were dissolved in UltraPure water and Mice for tissue isolation cDNA synthesis was performed using QuantiTect Re- Proximal small intestine was isolated from C57BL/6 verse Transcription Kit (Qiagen). qPCR reactions were mice of both sexes, aged between 3 and 6 months in all performed using TaqMan Universal Master Mix II (no experiments. UNG), pre-designed TaqMan probes (Additional file 9: Table S5), and 50 ng of sample cDNA (LifeTech). Reac- Bacteria strains tions were performed using an Applied Biosystems Cells were stored at –80 ºC and grown as follows. E. coli 7900HT system. qPCR results were analyzed using RQ strain MG1655 was grown overnight in LB. For experi- manager 1.2 software to obtain CT values used for rela- ments, overnight cultures of MG1655 were resuspended tive quantification to the housekeeping gene Hprt. in M9 supplemented with 0.4% glucose and 0.2% CASa- mino acids. L. lactis strain MG1363 was grown in M17 Confocal imaging of whole cell clusters media supplemented with 0.5% glucose, and E. faecalis ISC-enriched cell clusters (ENR + CV) suspended in 40 μL strain V583 was grown in Brain Heart Infusion media. of Matrigel were seeded onto round coverslips inside a 24-well plate. Cells were treated with ENR + CD, ENR + Crypt culture, enrichment, and differentiation CV, or ENR as previously described. At day 6, organoids Small intestinal crypts were cultured as previously de- were rinsed (PBS0 3X) and fixed on the coverslips by incu- scribed [59]. Briefly, crypts were resuspended in basal bating with 4% paraformaldehyde (PFA) for 30 min at culture medium (Advanced DMEM/F12 with 2 mM room temperature (RT). Gels were blocked and perme- GlutaMAX and 10 mM HEPES; Thermo Fisher Scien- abilized by incubating at RT for 1 hour with 0.1% Triton tific) at a 1:1 ratio with Corning™ Matrigel™ Membrane X-100 and 5% Powerblock in PBS0. Organoids were Matrix – GFR (Fisher Scientific) and plated at the center stained for DEFA and LYZ by incubating with rat of each well of 24-well plates. Following Matrigel anti-mouse Crp1 (Ayabe Lab clone 77-R63, 5 μg/mL, 50X) polymerization, 500 μL of small intestinal crypt culture and rabbit anti-human Lyz (Dako, RRID: AB_2341230, medium (basal media plus 100X N2 supplement, 50X B27 200X) primary antibodies diluted to 10 μg/mL in staining supplement; Life Technologies, 500X N-acetyl-L-cysteine; solution (0.1% Triton X-100 and 10X Powerblock in PBS0) Mead et al. BMC Biology (2018) 16:62 Page 17 of 24 overnight at 4 °C, followed by secondary antibodies Alexa for 20 min to dissociate into single cells. Dissociated Fluor 647 anti-Rabbit IgG (RRID: AB_2563202, 400X) and cells were centrifuged at 300 g for 3 min at 4 °C. The Alexa Fluor 488 anti-Rat IgG (RRID: AB_2563120, 400X) pellet was resuspended in FACS buffer (1% FBS in PBS, diluted in staining solution for 1 h at RT. Actin was stained Thermo Fisher Scientific) and strained into a 5 mL filter with Alexa Fluor 555 Phalloidin (40X) for 20 min, followed cap tube using a 40 μm filter. The cell suspension was by staining of the nucleus with 3 μMDAPIfor 5min. transferred to a flow prep microcentrifuge tube and Coverslips were mounted onto slides with Vectashield and centrifuged at 300 rcf for 3 min. Cell pellets were resus- imaged within 5 days using an Olympus FV2000 confocal pended in a Zombie violet dye (BioLegend, 100X) in microscope. Whole organoid confocal microscopy images FACS buffer for viability staining followed by 1% PFA were processed and analyzed using ImageJ. To determine fixation for 20 min at RT. Pellets were permeabilized thePC puritypercentage, theImageJPoint Picker plugin for 20 min at RT with staining buffer (0.5% Tween-20 wasused to count the numberofnucleito determine the in FACS buffer, Sigma), and co-stained with rabbit total number of cells and to count the number of DEFA- anti-human FITC-Lyz (Dako, RRID: AB_578661, 100X) and LYZ-containing PCs across all z-slices. To investigate and rat anti-mouse APC-CD24 (Biolegend, RRID: cell polarity in whole organoids, individual cells were AB_2565650, 100X) antibodies diluted in staining selected using ImageJ and mean area intensity within buffer for 45 min at RT. Flow cytometry was per- selected cell areas was computed in each z-slice through- formed using a BD LSR II HTS (BD; Koch Institute out the depth of the image across every channel imaged. Flow Cytometry Core at MIT). Initial settings and laser voltages were determined with unstained, single chan- High-resolution single-cell imaging nel stains or secondary-only controls (data not shown). Cell clusters were harvested and rinsed (basal culture Flow cytometry data was analyzed using FlowJo v10.7 media 3X) to remove Matrigel as previously described. software. Briefly, gating was performed as seen in Isolated clusters were resuspended in TrypLE Express Additional file 3: Figure S1F by removing doubles and and incubated at 37 °C for 20 min to dissociate into debris, then selecting the BV421 (viable) cell popula- single cells, then rinsed (basal culture media 2X) and tion; within this population, gating was based on LYZ resuspended in PBS containing magnesium and calcium. and CD24 populations. Pre-coated poly-L-lysine coverslips (Fisher Scientific) were placed into wells of a 24-well plate, a cell suspen- Lysozyme functional secretion assay sion containing approximately 50,000 cells per well was Lysozyme secretion was measured using a Lysozyme added to each well, and the plate was centrifuged at 700 Assay Kit (EnzChek; Thermo Fisher). Briefly, cells sus- rcf for 5 min. PBS supernatant was removed from the pended in Matrigel in 24-well plates were washed (basal wells, and the cells attached to the coverslips were fixed culture media 3X) and either supplemented with 500 μL by incubating with 4% PFA for 30 min at RT. After each of basal culture media or basal culture media plus step, cells were rinsed (PBS 2–5 min 3X). Cells were 10 μM CCh (Sigma Aldrich) for 24 h at 37 °C. Following blocked and permeabilized by incubating at RT for stimulation, culture plates were spun at high speed (> 30 min with permeabilization solution and stained with 2000 g) for 5 min at RT to pellet cell debris and loose DEFA and LYZ by incubating with rat anti-mouse Crp1 Matrigel, and 25 μL of conditioned supernatant was (Ayabe Lab clone 77-R63, 5 μg/mL, 50X) and rabbit removed from the top of each well and quantified as per anti-human Lyz (Dako, RRID: AB_2341230, 200X) pri- the manufacturer’s protocol. mary antibodies diluted in staining solution overnight at 4 °C. Secondary antibodies Alexa Fluor 647 anti-Rabbit Quantification of cell viability, apoptosis, and cytotoxicity IgG (RRID: AB_2563202, 400X) and Alexa Fluor 488 To track proliferation and cell viability, DNA content anti-Rat IgG (RRID: AB_2563120, 400X) diluted in stain- was quantified over the course of differentiation and ing solution were incubated with the coverslips for 1 h CCh stimulation using a CyQUANT Cell Proliferation at RT. Actin was stained with Alexa Fluor 555 Phalloidin Assay Kit (Thermo Fisher) as per the manufacturer’s incubated for 20 min at RT, and the nucleus was stained protocol. Briefly, culture media was aspirated from each with DAPI by incubating at RT for 5 mins. Coverslips well, and the wells were washed (PBS 3X). Gels were were mounted on to slides with Vectashield and imaged then mechanically dissociated into PBS, the contents within 48 h using an Applied Precision DeltaVision transferred into a Falcon tube, centrifuged at 300 rcf for Microscope. 3 min at 4 °C, and the pellet resuspended in PBS to wash. Tubes were centrifuged at 300 rcf for 5 min at Flow cytometry 4 °C, and the pellet was resuspended in 1 mL of assay Cell clusters were isolated from Matrigel as previously working solution (20X cell-lysis buffer, 400X GR dye in described and resuspended in TrypLE Express at 37 °C DI water); 200 μL of samples and DNA standards were Mead et al. BMC Biology (2018) 16:62 Page 18 of 24 plated in triplicate in a black 96-well plate, shaken for concentration of 1.0% and subsequently desalted on 5 min, then fluorescence was measured on a plate 10 mg OASIS HLB solid phase columns (Waters). reader (480 nm/520 nm). From each condition (n = 8), 50 μg aliquots of the Ng For ENR/ENR + CD co-culture, ISC-enriched orga- KD dried tryptic peptides were reconstituted in 100 mM noids (ENR + CV) were differentiated in ENR and ENR HEPES at pH 8.0 to a final concentration of 1.0 mg/mL. + CD and isolated as previously described. The cell The peptides were labeled with TMT-10 isobaric mass pellets were counted and resuspended in basal culture tag reagent according to the manufacturer’s instructions medium, mixed at 0:100, 25:75, 50:50, 75:25, and 100:0% (ThermoFisher Scientific). The peptides were labeled at ENR:ENR + CD ratios (number of clusters), and plated a 1:8 ratio of peptide to TMT reagent, followed by 1 h as previously described in Matrigel in a 96-well plate at incubation at RT with bench-top shaking at 850 rpm. approximately 50 clusters per well in ENR media. After After incubation, a 1.0 μg aliquot of labeled tryptic 24 and 48 h of co-culture, viability versus cytotoxicity peptide was removed from each labeled condition, and caspase activation were assessed using ApoTox-Glo desalted with C18 stage tips, and analyzed via LC-MS/ Triplex Assay (Promega) according to the manufac- MS using a Thermo Fisher Q Exactive Plus Hybrid Mass turer’s protocol. Briefly, 20 μLof ‘V/C reagent’ (10 μL Spectrometer (QE-Plus) coupled to a Thermo Fisher each of GF-AFC and bis-AAF-R110 substrates in 2.0 mL EASY-nLC 1000 liquid chromatograph to ensure iso- of assay buffer) were added to all wells and mixed by baric label incorporation ≥ 95%. An additional 1.0 μgof orbital shaking at 500 rpm for 30 s. After 30 min of labeled tryptic peptide was removed from each channel, incubation at 37 °C, fluorescence was measured on a mixed together, desalted on a C18 stage tip, and ana- plate reader (400 nm/505 nm for viability and 485 nm/ lyzed via LC-MS to ensure equal relative protein loads. 520 nm for cytotoxicity). Caspase-Glo 3/7 reagent (100 μL) During these quality control steps, the labeled peptides was then added to all wells and mixed by orbital shaking at were stored, unquenched, at –80 °C. After validation, 500 rpm for 30 s. After 30 min of incubation at RT, lumi- each channel was quenched with a 5% hydroxylamine nescence was measured on a plate reader. solution to a final sample concentration of 0.3% to quench any unbound isobaric tags. The corresponding eight channels were mixed together for a total amount Bacteria co-culture of 400 μg of labeled tryptic peptides. The labeled peptide For bacteria co-culture, ISC-enriched cells (ENR + CV) mixture was dried down in a speedvac and subsequently were differentiated in ENR and ENR + CD as previously desalted on 30 mg of OASIS HLB solid phase column described. After 6 days of differentiation, cell clusters (Waters). were isolated as previously described. The cell pellet was The dried, labeled peptides were fractionated into 24 resuspended in basal culture medium and plated in sus- fractions by basic reversed-phase using an Agilent pension in a 96-well plate at approximately 150 clusters Zorbax 300 Å, 4.6 mm × 250 mm Extend-C18 column per well. A 1:1 volume of bacteria in respective media on an Agilent 1100 Series HPLC instrument (Agilent (see ‘Bacterial strains’, above; in exponential growth, as Technologies) to decrease sample complexity and in- confirmed by plate reader OD) was added, and bacterial crease the dynamic range of detection. Solvent A (2% growth was measured by serial plating (CFU) after a 4 h acetonitrile, 5 mM ammonium formate, pH 10), and a incubation. Results for bacteria co-culture were normal- non-linear increasing concentration of solvent B (90% ized to no cell (bacteria only) controls. acetonitrile, 5 mM ammonium formate, pH 10) was used as the mobile phase with a flow rate of 1 mL/min Mass spectrometry proteomics sample preparation, through the column. A non-linear gradient with increas- sequencing, and quantification ing percentages of solvent B with four different slopes Organoid cell pellets were isolated from Matrigel with was used (0% for 7 min; 0% to 16% in 6 min; 16% to mechanical dissociation and washed (cold PBS 5X) to 40% in 60 min; 40% to 44% in 4 min; 44% to 60% in remove residual extracellular protein. Proteins were ex- 5 min; 60% for 14 min), and the eluted peptides were tracted from cell pellets with 8 M urea (Sigma), reduced collected in a Whatman polypropylene 2 mL 96-well with 5 mM DTT (Thermo Fisher Pierce) for 45 min, plate (Whatman). The 96 fractions were concatenated alkylated with 10 mM IAA (Sigma) for 45 min in the down to 25 fractions. dark, and double digested with both Lysyl Endopeptidase The global proteome (25 fractions) was analyzed by ‘LysC’ (Wako) and trypsin (Promega) overnight at RT. A LC-MS/MS using the same system described above. –1 small aliquot of cellular lysate was removed from each Peptides were separated at a flow rate of 200 nL min sample for protein quantification via the Pierce™ BCA on a capillary column (Picofrit with a 10-μm tip opening Protein Assay Kit (Pierce). After proteolytic digestion, and 75 μm diameter, New Objective, PF360-75-10-N-5) the samples were quenched using formic acid to a final packed at the Broad Institute with 20 cm of C18 1.9 μm Mead et al. BMC Biology (2018) 16:62 Page 19 of 24 silica beads (1.9 μm ReproSil-Pur C18-AQ medium, Dr. Therefore, in this step of normalization, we normalized Maisch GmbH, r119.aq). Injected peptides were sepa- the median log2 ratio for each ratio column so that the –1 rated at a flow rate of 200 nL min with a linear 84 min median log2 ratio was zero. To robustly and confidently gradient from 100% solvent A (3% acetonitrile, 0.1% for- detect real differential peptides and proteins in the mic acid) to 30% solvent B (90% acetonitrile, 0.1% formic TMT-labeled experiment, we performed a moderated t acid), followed by a linear 9 min gradient from 30% solv- test [63, 64]. Unlike the standard t test, which is not ent A to 90% solvent B for a total of 110 min. The robust for small numbers of samples, the moderated t test QE-Plus instrument was operated in the data-dependent uses an empirical Bayes approach that ‘moderates’ vari- mode acquiring higher-energy collisional dissociation ance estimates for peptides (i.e., shrunk towards a com- tandem mass spectrometry (HCD MS/MS) scans (Reso- mon value), thereby significantly improving the stability of lution = 35,000) for TMT-10 on the 12 most abundant variance estimates for individual peptides. The p values re- ions using an MS1 ion target of 3 × 10 ions and an MS2 ported by the moderated t test were adjusted for multiple target of 5 × 10 ions. The maximum ion time used for testing using the Benjamini–Hochberg FDR method [64]. the MS/MS scans was 120 ms; the HCD-normalized collision energy was set to 31; the dynamic exclusion Proteome pathway and network analysis time was set to 20 s, and the peptide-match preferred Using the identified and quantified proteins from the setting was enabled. TMT-10 labeling experiment, multiple pathway and net- work analyses were performed. Sample correlations were Protein and peptide identification and quantification represented as r-values and determined using GraphPad Peptide spectrum matching and protein identification was Prism version 7.0a. To elucidate potential transcriptional performed using Agilent Technologies SM software pack- drivers of proteome structure, we performed GSEA age (developed at the Broad Institute). In SM, FDRs are v3.0b2 [37, 38] using the full rank-ordered proteome calculated at three different levels, namely spectrum, dis- against the transcription factor target gene set database tinct peptide, and distinct protein. Peptide FDRs are calcu- (v5.2 MSigDB) [39], then performed enrichment map lated in SM using essentially the same pseudo-reversal visualization using GSEA-P-based implementation and strategy evaluated by Elias and Gygi [60] and shown to Cytoscape v3.4.0 [40, 41], with a moderately conservative perform the same as library concatenation. A false distinct cutoff (p < 0.005 and FDR < 0.075) and an overlap coeffi- protein ID occurs when all the distinct peptides that group cient of 0.2. To assess the functional and compartmen- together to constitute a distinct protein have a deltaFor- tal functions associated with the ENR + CD-enriched wardReverseScore ≤ 0. We adjusted the settings to provide proteome and ENR-enriched proteome, we used DA- peptide FDR of 1–2% and protein FDR of 0–1%. SM also VID v6.8 [65, 66] and the gene ontology database, look- carries out sophisticated protein grouping using the ing only at experimentally verified associations within methods previously described [61]. Only proteins with biological processes, cellular compartments, and mo- more than two peptides and at least two TMT ratios in lecular function against a background set of all 8015 each replicate were counted as being identified and quan- quantified proteins. tified. Additionally, we added the capability to flag poten- tially unreliable TMT quantification results based on Single-cell RNA-sequencing detection of more than one precursor in the selection win- A single-cell suspension was obtained from organoids cul- dow for MS/MS. The precursor ion flagging was similar tured under ENR + CV, ENR, and ENR + CD conditions to that recently reported [62], but was carried out post for 6 days as described above. We utilized the Seq-Well data acquisition. As an output, SM generates protein and platform for massively parallel scRNA-seq to capture tran- peptide reports for downstream differential regulation, scriptomes of single cells on barcoded mRNA capture pathway, and network analysis. Prior to comprehensive beads. Full methods on implementation of this platform differential marker, pathway, and network analysis with are available in Gierahn et al. [28]. In brief, 20,000 cells the SM-generated protein reports, we ensured that the from one organoid condition were loaded onto one array data was of high quality and has been properly normal- containing 100,000 barcoded mRNA capture beads. The ized. The first level of normalization was accomplished by loaded arrays containing cells and beads were then sealed guaranteeing that an equivalent amount of peptide (50 μg using a polycarbonate membrane with a pore size of per) was labeled for each of the 10 TMT channels. Once 0.01 μm, which allows for exchange of buffers but retains the SM reports were generated, we calculated the median biological molecules confined within each microwell. ratios for each of the channels where the denominator of Subsequent exchange of buffers allows for cell lysis, tran- the ratio was a predetermined TMT channel signifying the script hybridization, and bead recovery before performing control condition. The underlying assumption was that reverse transcription en masse. Following reverse tran- the null distribution is centered at zero in log2 space. scription and exonuclease treatment to remove excess Mead et al. BMC Biology (2018) 16:62 Page 20 of 24 primers, PCR amplification was carried out using KAPA total number of ENR + CV, ENR, and ENR + CD cells in- HiFi PCR Mastermix with 2000 beads per 50 μLreaction cluded in the analysis was 985, 2544, and 2382, respect- volume. Six libraries (totaling 12,000 beads) were then ively, with quality metrics for nGene, nUMI, and pooled and purified using Agencourt AMPure XP beads percentage of ribosomal and mitochondrial genes re- (Beckman Coulter, A63881) by a 0.6X SPRI followed by a ported in Additional file 8: Figure S4. We then per- 0.7X SPRI and quantified using Qubit hsDNA Assay formed principal component analysis over the list of (Thermo Fisher). Libraries were constructed using the variable genes. For both clustering and t-stochastic Nextera Tagmentation method on a total of 800 pg of neighbor embedding (tSNE), we utilized the first 12 pooled cDNA library from 12,000 recovered beads. Tag- principal components based on the elbow method, as mented and amplified sequences were purified at a 0.6X upon visual inspection of genes contained within, each SPRI ratio yielding library sizes with an average distribu- contributing to important biological processes of intes- tion of 650–750 base pairs in length as determined using tinal cells. We used FindClusters with a resolution of the Agilent hsD1000 Screen Tape System (Agilent Gen- 1.35 and 1000 iterations of tSNE to identify 14 clusters omics). Arrays were sequenced with an Illumina 75 Cycle across the three input samples. To identify genes NextSeq500/550v2 kit at a final concentration of 2.8 pM. which defined each cluster, we performed a ROC test The read structure was paired end with Read 1 starting implemented in Seurat with a threshold set to an from a custom read 1 primer containing 20 bases with AUC of 0.60. a 12 bp cell barcode and 8 bp UMI and Read 2 being 50 bases containing transcript information (Additional Transcriptional scoring file 10: Table S6 for primers used). To determine the fractional contribution to a cell’s tran- scriptome of a gene list, we summed the total log(scaled Single-cell RNA-sequencing computational pipelines and UMI + 1) expression values for genes within a list of analysis interest and divided by the total amount of scaled UMI Read alignment was performed as in Macosko et al. [67]. detected in that cell giving a proportion of a cell’s tran- Briefly, for each NextSeq sequencing run, raw sequen- scriptome dedicated to producing those genes. From our cing data was converted to demultiplexed FASTQ files proteomic screen, we took a list of upregulated proteins using bcl2fastq2 based on Nextera N700 indices corre- (249) or downregulated proteins (212) that were de- sponding to individual samples/arrays. Reads were then tected within our single-cell RNA-sequencing data. To aligned to mm10 genome using the Galaxy portal main- determine the relationship to in vivo PCs and EECs, we tained by the Broad Institute for Drop-Seq alignment took reference data from two Seq-Well experiments run using standard settings. Individual reads were tagged on epithelial cells dissociated from the ileal region of the according to the 12 bp barcode sequencing and the 8 bp small intestine of two C57BL/6 J mice run in separate UMI contained in Read 1 of each fragment. Following experiments. The ileum was first rinsed in 30 mL of ice alignment, reads were binned onto 12 bp cell barcodes cold PBS and allowed to settle. The segment was then and collapsed by their 8 bp UMI. Digital gene expression sliced with scissors and transferred to 10 mL of epithe- matrices (e.g., cell by gene tables) for each sample were lial cell solution (HBSS Ca/Mg-Free 10 mM EDTA, obtained from quality filtered and mapped reads and 100 U/mL penicillin, 100 μg/mL streptomycin, 10 mM UMI-collapsed data, were deposited in GSE100274, and HEPES, 2% FCS (ThermoFisher)) freshly supplemented were utilized as input into Seurat (https://github.com/ with 200 μL of 0.5 M EDTA. The epithelial separation satijalab/seurat) for further analysis [68]. from the underlying lamina propria was performed for To analyze ENR + CV, ENR, and ENR + CD organoids 15 min at 37 °C in a rotisserie rack with end-over-end together, we merged UMI matrices across all genes de- rotation. The tube was then removed and placed on ice tected in any condition and generated a matrix retaining immediately for 10 min before shaking vigorously 15 all cells with at least 1000 UMI detected. This table was times. Visual macroscopic inspection of the tube at this then utilized to setup the Seurat object in which any cell point yielded visible epithelial sheets, and microscopic with at least 400 unique genes was retained and any examination confirmed the presence of single-layer gene expressed in at least five cells was retained. The sheets and crypt-villus structures. The epithelial fraction object was initiated with log-normalization, scaling, and was spun down at 400 g for 7 min and resuspended in centering set to True. Before performing dimensionality 1 mL of epithelial cell solution before transferring to a reduction, data was subset to include cells with less than 1.5 mL Eppendorf tube to minimize the time spent 8000 UMI, and a list of 1676 most variable genes was centrifuging. Cells were spun down at 800 g for 2 min generated by including genes with an average normal- and resuspended in TrypLE Express for 5 min in a 37 °C ized and scaled expression value greater than 0.14 and bath followed by gentle trituration with a P1000 pipette. with a dispersion (variance/mean) greater than 0.4. The Cells were spun down at 800 g for 2 min and Mead et al. BMC Biology (2018) 16:62 Page 21 of 24 resuspended in ACK lysis buffer (ThermoFisher) for Additional file 7: Figure S3 and Additional file 8:FigureS4; 3 min on ice to remove red blood cells and dying cells. n = 8 single-well replicates from one and five biological Cells were spun down at 800 g for 2 min and resus- donors for data reported in Fig. 5a, b and c, respectively; pended in 1 mL of epithelial cell solution and placed on n = 13 co-culture well replicates randomly selected with- ice for 3 min before triturating with a P1000 pipette and out replacement from four donors for data reported in filtering into a new Eppendorf through a 40 μmcell Fig. 5d; n = 6 well replicates (two per three biological strainer (Falcon/VWR). Cells were spun down at 800 g for donors) in Fig. 5e; n = 3 biological donors in Fig. 6c,and 2 min and then resupended in 200 μL of epithelial cell so- n = 5 biological donors in Fig. 6d. For scores in single-cell lution and placed on ice for counting. Single-cell RNA-seq data, we report effect sizes in addition to statistical signifi- data was then generated as described in the ‘Single-cell cance as an additional metric for the magnitude of the RNA-sequencing’ and ‘Single-cell RNA-sequencing com- effect observed. The calculation was performed as Cohen’s putational pipelines and analysis’ sections in Methods.To d, where effect size d =(Mean -Mean )/(SD pooled). 1 2 generate PC and EEC signatures, we ran unbiased SNN-graph based clustering, performed a ROC test, identi- Additional files fied the two mature PC and EEC clusters, and reported all genes with an AUC above 0.60, using all genes with an Additional file 1: Table S1. Derived gene list of the top defining genes from in vivo ileal small intestine PCs and EECs captured on the Seq-Well AUC above 0.65 for scoring within each cluster (gene lists platform. (XLSX 355 kb) in Additional file 1: Table S1), representing any gene with Additional file 2: Table S2. Reference gene lists used in single-cell enrichment in PCs and EECs. These lists capture genes analyses. (XLSX 24 kb) which are enriched in PC (Lyz-high) and EECs (Chga-high) Additional file 3: Figure S1. Image analysis of cell clusters and flow and separate them from the rest of the cells present in in- cytometry. A Bright-field microscopy after 6 days of ENR + CD culture shows annular morphology and darkened lumen of cell clusters consistent with testinal epithelium. For pathway analysis, we inspected cu- presence of granule-rich cells. B Percentage of total cells that are LYZ+ and rated gene lists deposited in the GSEA platform and used DEFA+ following 6 days of ENR, ENR + CV, and ENR + CD culture (from cell KEGG-derived Wnt and Reactome-derived Notch and Re- counting of whole clusters) (n = 3 minimum biological replicates, one-way ANOVA with multiple comparison test versus ENR, **** adj. p < 0.0001). spiratory Electron Transport Chain signatures (Additional C Collapsed z-stack of whole cluster with individual cells highlighted (1–3) file 2: Table S2). In vivo transcription factors for PCs and following 6 days of ENR + CD, stained for LYZ and DEFA and counterstained EECs were determined by matching the PC and EEC signa- with DAPI and for actin (phalloidin). (1–3) Normalized mean-area intensity versus z-axis depth profiles of representative individual LYZ+/DEFA+ ture gene sets with transcription factors from the Riken co-staining cells. D Representative flow cytometry of ENR and ENR + CD Transcription Factor Database (TFdb - http://genome.g- at 6 days with distinct populations of CD24+ and LYZ+ cells indicative sc.riken.jp/TFdb/), and then including only those TFs of phenotypic PCs. E Representative gating for flow cytometry, including removal of doublets and non-viable cells in final gating. F Percentage of which were robustly identified in the proteome dataset. viable cells (membrane impermeable) over time of ENR versus ENR + CD culture. (PDF 3217 kb) Quantification and statistical analysis Additional file 4: Figure S2. Proteomic pipeline, sample-to-sample All statistical analyses were performed using GraphPad comparison, and insights from the in vitro PC proteome. A Schematic of proteomic analysis for samples: culture, collection, lysis, reduction and Prism v7.0a, Seurat implemented in RStudio, and Agilent alkylation, proteolytic digestion, labeling of peptides with isobaric mass tag Technologies Spectrum Mill software package. All graphs reagents (Tandem Mass Tags, TMT10-plex; Thermo), off-line fractionation by with n > 6 show mean ± SEM, unless otherwise noted, basic reverse phase chromatography, analysis of fractions by LC-MS/MS, identification of peptides and proteins using Spectrum Mill software whereas graphs with n < 6 show mean and individual repli- (Agilent), and statistical analysis of the resulting data (moderated t test) to cate values. Unpaired two-tail t test and two-way ANOVA identify confidently differential proteins. B Proteome sample correlation with Dunnett’s multiple comparison test (reported as adj. p between all biological (n = 2) and technical (n = 2/biological) replicates. C ENR + CD-enriched proteins are well-annotated in the gene ontology value) were used to assess statistical significance as appro- (GO) database and show robust enrichment for functions and compartments priate and, unless otherwise noted, * indicates p <0.05, ** p of secretory cells determined by fold enrichment vs. FDR using DAVID. D ENR- < 0.01 *** p < 0.001, **** p < 0.0001, and ns non-significant. enriched proteins are well annotated in the GO database and show enrichment for functions and compartments of transcriptionally and In each experiment, tissues were isolated from multiple translationally active cells determined by fold enrichment vs. FDR mice housed in the same facility with each mouse providing using DAVID. e GSEA enrichment map of transcription factors linked tissue designated as a distinct biological donor: n =3 to ENR + CD- and ENR-enriched proteins following a moderately conservative cutoff of p < 0.005, FDR < 0.075, and overlap coefficient donor-averaged values of four technical replicates for data of 0.2. (PDF 483 kb) reported in Fig. 2b; n = 3 donor of two technical replicates Additional file 5: Table S3. Detected and quantified in vitro Proteome. for data reported in Fig. 2e (with exception of day 8, n =2 (XLSX 919 kb) donor-averaged values of two technical replicates) and Additional file 6: Table S4. Complete list of DAVID enrichments. Additional file 3:FigureS1F; n = 4 (two technical replicates (XLSX 51 kb) from two biological donors each) for data reported in Additional file 7: Figure S3. Quality metrics for single-cell RNA sequencing. A Total gene number of cells maintained in analyses with a lower cutoff of Figs. 2f–h, and Additional file 4:FigureS2; n = 1 biological n = 400 unique genes per cell. Total unique molecular identifiers (UMIs) used donor for in vitro data reported in Figs. 3 and 4 and Mead et al. BMC Biology (2018) 16:62 Page 22 of 24 Competing interests as the basis for cell-by-gene tables collapsed to UMI as input into Seurat with JMK, RSL, and XY hold equity in Frequency Therapeutics, a company that has the lower bound representing n = 400 unique genes and the upper bound an option to license IP generated by JMK, RSL, and XY and that may benefit representing 8000 UMIs. Note: Clusters ENR + CV-3, ENR + CV-4, and ENR-1 had financially if the IP is licensed and further validated. The interests of JMK, RSL, significantly higher levels of genes and UMIs and, intriguingly, were also the and XY were reviewed and are subject to a management plan overseen by three clusters with the highest levels of Lgr5 (see Fig. 5a), indicating that stem their institutions in accordance with their conflict of interest policies. cells may contain larger contents of RNA, as they are in a biosynthetic state before differentiation and maturation. B Violin plot of expression contribution to a cell’s transcriptome of mitochondrial and ribosomal Publisher’sNote genes across identified subsets. (PDF 1153 kb) Springer Nature remains neutral with regard to jurisdictional claims in Additional file 8: Figure S4. Signaling pathways and processes associated published maps and institutional affiliations. with in vitro PC enrichment (Additional file 9: Table S5 for reference gene lists). A Violin plot of expression contribution to a cell’s transcriptome of Wnt Author details pathway genes (Additional file 2: Table S2; activated by CHIR99021) across 1 Division of Engineering in Medicine, Department of Medicine, Brigham & clusters as percent of transcriptome. B Violin plot of expression contribution 2 Women’s Hospital, Harvard Medical School, Boston, MA, USA. Harvard-MIT to a cell’s transcriptome of Notch pathway genes (Additional file 2:Table S2; 3 Division of Health Sciences & Technology, Cambridge, MA, USA. Koch inhibited by DAPT) across clusters as percent of transcriptome. C Violin plot Institute for Integrative Cancer Research, MIT, Cambridge, MA, USA. Harvard of expression contribution to a cell’s transcriptome of respiratory electron 5 Stem Cell Institute, Cambridge, MA, USA. Broad Institute of Harvard and MIT, transport gene set (Additional file 2:Table S2) acrossclustersaspercent of 6 Cambridge, MA, USA. Institute for Medical Engineering and Science, MIT, transcriptome. (PDF 888 kb) 7 Cambridge, MA, USA. Department of Chemistry, MIT, Cambridge, MA, USA. 8 9 Additional file 9 : Table S5. TaqMan gene expression assays used for Department of Chemical Engineering, MIT, Cambridge, MA,, USA. Ragon qRT-PCR. (DOCX 13 kb) Institute of MGH, MIT and Harvard, Cambridge, MA, USA. Divisions of Infectious Diseases and Gastroenterology, Boston Children’s Hospital and Additional file 10: Table S6. SeqWell reverse transcription and library Harvard Medical School, Boston, MA, USA. Wyss Institute for Biologically preparation primers. (DOCX 13 kb) Inspired Engineering, Harvard University, Boston, MA, USA. Department of Biological Engineering, MIT, Cambridge, MA, USA. Synthetic Biology Center, MIT, Cambridge, MA, USA. Center for Microbiome Informatics and Acknowledgements Therapeutics, MIT, Cambridge, MA, USA. We thank the Koch Institute of MIT Core facilities for assisting with microscopy, flow cytometry, and pilot mass spectrometry experiments, and Received: 30 March 2018 Accepted: 8 May 2018 for the partial support from the Cancer Center Support (core) Grant P30- CA14051 from the NCI. We also thank R Blumberg, N Krupka, J Matute, and Y Shao for critically reviewing the manuscript. Furthermore, we appreciate the contribution of anti-crypitidin (DEFA) antibodies from Tokiyoshi Ayabe References (Hokkaido University). 1. Clevers H. Modeling development and disease with organoids. Cell. 2016; 165:1586–97. https://doi.org/10.1016/j.cell.2016.05.082 Funding 2. Prakadan SM, Shalek AK, Weitz DA. Scaling by shrinking: empowering This research was supported by both an Innovator award and Breakthrough single-cell “omics” with microfluidic devices. Nat Rev Genet. 2017;18:345–61. award to JMK from the Kenneth Rainin Foundation, the US National Institutes of 3. Haber AL, Biton M, Rogel N, Herbst RH, Shekhar K, Smillie C, et al. A single- Health (NIH) grant DE013023 to RL, NIH grant HL095722 to JMK, and a ‘shark tank’ cell survey of the small intestinal epithelium. Nature. 2017;551(7680):333–9. pilot grant from Brigham and Women’sHospitaltoJMK.BEM wassupported by 4. Grün D, Lyubimova A, Kester L, Wiebrands K, Basak O, Sasaki N, et al. Single- the Graduate Research Fellowship Program from the US National Science cell messenger RNA sequencing reveals rare intestinal cell types. Nature. Foundation (NSF). AKS was supported by the Searle Scholars Program, the 2015;525:251–5. Beckman Young. 5. The HCA Consortium. The human cell atlas white paper. 2017; https://www. Investigator Program, a Sloan Research Fellowship in Chemistry, NIHgrants humancellatlas.org/files/HCA_WhitePaper_18Oct2017.pdf. Accessed 29 May 2018. 1DP2OD020839, 2U19AI089992, 2R01HL095791, 1U54CA217377,2P01AI039671, 6. Tanay A, Regev A. Scaling single-cell genomics from phenomenology to 5U24AI118672, 2RM1HG006193, 1R33CA202820, 1R01HL126554,1R01DA046277 mechanism. Nature. 2017;541:331–8. and 1R01AI138546, and Bill and Melinda Gates Foundationgrants OPP1139972, 7. Satija R, Shalek AK. Heterogeneity in immune responses: From populations OPP1137006, and OPP1116944. JOM is a Damon Runyon Fellow supported by to single cells. Trends Immunol. 2014;35:219–29. https://doi.org/10.1016/j.it. the HHMI and Damon Runyon Cancer Research Foundation DRG-2274-16. 2014.03.004 8. Lancaster MA, Knoblich JA. Organogenesis in a dish: Modeling development and disease using organoid technologies. Science. 2014;345:1247125. Availability of data and materials 9. Schwank G, Koo BK, Sasselli V, Dekkers JF, Heo I, Demircan T, et al. The accession number for the proteomics data reported in this paper is Functional repair of CFTR by CRISPR/Cas9 in intestinal stem cell organoids MassIVE: MSV000081369. The accession number for the single-cell RNA of cystic fibrosis patients. Cell Stem Cell. 2013;13:653–8. https://doi.org/10. sequencing data reported in this paper is GEO: GSE100274. 1016/j.stem.2013.11.002 10. Drost J, van Boxtel R, Blokzijl F, Mizutani T, Sasaki N, Sasselli V, et al. Use of Authors’ contributions CRISPR-modified human stem cell organoids to study the origin of Conceptualization, BEM, JOM, XY, AKS and JMK; Methodology, BEM, XY, JMK, mutational signatures in cancer. Science. 2017;238:eaao3130. JOM, AKS, RA; Formal Analysis, BEM, JOM, MJS, MAM, Investigation, BEM, APB, 11. Molodecky NA, Soon IS, Rabi DM, Ghali WA, Ferris M, Chernoff G, et al. LEL, JOM, MJS, DAA, PB, TKH, MHW, SRN; Resources, SRN, RL, JMK, SAC, AKS, Increasing incidence and prevalence of the inflammatory bowel diseases Writing – Original Draft, BEM, LEL, JOM, APB; Writing – Review & Editing, with time, based on systematic review. Gastroenterology. 2012;142:46–54. BEM, JMK, LEL, JOM, PB, XY, RA, SAC, JJC, AKS, RL; Visualization, BEM, JOM, e42; quiz e30 LEL, APB; Supervision, JMK, BEM, JJC, AKS, SAC; Project Administration, JMK, 12. Wehkamp J, Salzman NH, Porter E, Nuding S, Weichenthal M, Petras RE, BEM; Funding Acquisition, JMK, AKS, RL, BEM. All authors read and approved et al. Reduced Paneth cell α-defensins in ileal Crohn’s disease. Proc Natl the final manuscript. Acad Sci U S A. 2005;102:18129–34. 13. Ireland H, Houghton C, Howard L, Winton DJ. Cellular inheritance of a Ethics approval and consent to participate Cre-activated reporter gene to determine Paneth cell longevity in the Isolation of tissues, housing, and maintenance of animal colonies were all murine small intestine. Dev Dyn. 2005;233:1332–6. performed under the approval and guidelines of the MIT Committee on 14. Sato T, van Es JH, Snippert HJ, Stange DE, Vries RG, van den Born M, et al. Animal Care or the Harvard Digestive Disease Center and Brigham and Paneth cells constitute the niche for Lgr5 stem cells in intestinal crypts. Women’s Gnotobiotic Core. Nature. 2011;469:415–8. Mead et al. BMC Biology (2018) 16:62 Page 23 of 24 15. Clevers HC, Bevins CL. Paneth cells: maestros of the small intestinal crypts. 40. Merico D, Isserlin R, Stueker O, Emili A, Bader GD. Enrichment map: a Annu Rev Physiol. 2013;75:289–311. network-based method for gene-set enrichment visualization and 16. Xavier RJ, Podolsky DK. Unravelling the pathogenesis of inflammatory bowel interpretation. PLoS One [Internet]. 2010;5:e13984. Available from: disease. Nature. 2007;448:427–34. http://www.ncbi.nlm.nih.gov/pubmed/21085593 17. Khor B, Gardet A, Xavier RJ. Genetics and pathogenesis of inflammatory 41. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. bowel disease. Nature. 2011;474:307–17. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504. 18. Liu T-C, Gurram B, Baldridge MT, Head R, Lam V, Luo C, et al. Paneth cell 42. Stringari C, Edwards RA, Pate KT, Waterman ML, Donovan PJ, Gratton E. defects in Crohn’s disease patients promote dysbiosis. JCI Insight. 2016;1:1–15. 19. Adolph TE, Tomczak MF, Niederreiter L, Ko H-J, Böck J, Martinez-Naves E, Metabolic trajectory of cellular differentiation in small intestine by Phasor et al. Paneth cells as a site of origin for intestinal inflammation. Nature. Fluorescence Lifetime Microscopy of NADH. Sci Rep. 2012;2:568. 2013;503:272–6. 43. Rodríguez-Colman MJ, Schewe M, Meerlo M, Stigter E, Gerrits J, Pras-Raves 20. Kobayashi KS, Chamaillard M, Ogura Y, Henegariu O, Inohara N, Nuñez G, M, et al. Interplay between metabolic identities in the intestinal crypt et al. Nod2-dependent regulation of innate and adaptive immunity in the supports stem cell function. Nature. 2017;543(7645):424–7. intestinal tract. Science. 2005;307:731–4. 44. Ayabe T, Satchell DP, Wilson CL, Parks WC, Selsted ME, Ouellette AJ. 21. Kaser A, Blumberg RS. ATG16L1 Crohn’s disease risk stresses the Secretion of microbicidal alpha-defensins by intestinal Paneth cells in endoplasmic reticulum of Paneth cells. Gut. 2014;63(7):1038–9. response to bacteria. Nat Immunol. 2000;1:113–8. 45. Kanamori M, Konno H, Osato N, Kawai J, Hayashizaki Y, Suzuki H. A genome- 22. Kaser A, Lee A-H, Franke A, Glickman JN, Zeissig S, Tilg H, et al. XBP1 links wide and nonredundant mouse transcription factor database. Biochem ER stress to intestinal inflammation and confers genetic risk for human Biophys Res Commun. 2004;322:787–93. inflammatory bowel disease. Cell. 2008;134:743–56. 23. Farin HF, Karthaus WR, Kujala P, Rakhshandehroo M, Schwank G, Vries RGJ, 46. Jia SN, Lin C, Chen DF, Li AQ, Dai L, Zhang L, et al. The transcription factor et al. Paneth cell extrusion and release of antimicrobial products is directly p8 regulates Autophagy in response to palmitic acid stress via a controlled by immune cell-derived IFN-γ. J Exp Med. 2014;211:1393–405. mammalian target of rapamycin (mTOR)-independent signaling pathway. J 24. Wilson SS, Tocchi A, Holly MK, Parks WC, Smith JG. A small intestinal Biol Chem. 2016;291:4462–72. organoid model of non-invasive enteric pathogen-epithelial cell 47. Grasso D, Bintz J, Lomberk G, Molejon MI, Loncle C, Garcia MN, et al. Pivotal interactions. Mucosal Immunol. 2014;8:1–10. role of the chromatin protein Nupr1 in Kras-induced senescence and 25. Foulke-Abel J, In J, Yin J, Zachos NC, Kovbasnjuk O, Estes MK, et al. Human transformation. Sci Rep. 2015;5:17549. enteroids as a model of upper small intestinal ion transport physiology and 48. Cano CE, Hamidi T, Sandi MJ, Iovanna JL. Nupr1: The Swiss-knife of cancer. pathophysiology. Gastroenterology. 2016;150:638–49. e8 J Cell Physiol. 2011;226:1439–43. 49. Imielinski M, Baldassano RN, Griffiths A, Russell RK, Annese V, Dubinsky M, 26. Moon C, VanDussen KL, Miyoshi H, Stappenbeck TS. Development of a primary et al. Common variants at five new loci associated with early-onset mouse intestinal epithelial cell monolayer culture system to evaluate factors inflammatory bowel disease. Nat Genet. 2009;41:1335–40. that modulate IgA transcytosis. Mucosal Immunol. 2013;7:818–28. 27. Basak O, Beumer J, Wiebrands K, Seno H, van Oudenaarden A, Clevers H. 50. Santofimia-Castaño P, Rizzuti B, Pey ÁL, Soubeyran P, Vidal M, Urrutia R, Induced quiescence of Lgr5+ stem cells in intestinal organoids enables et al. Intrinsically disordered chromatin protein NUPR1 binds to the differentiation of hormone-producing enteroendocrine cells. Cell Stem Cell. C-terminal region of Polycomb RING1B. Proc Natl Acad Sci. 2017; 2017;20:177–90. e4 https://doi.org/10.1073/pnas.1619932114. 28. Gierahn TM, Wadsworth MH, Hughes TK, Bryson BD, Butler A, Satija R, et al. 51. Neira JL, Bintz J, Arruebo M, Rizzuti B, Bonacci T, Vega S, et al. Identification Seq-Well: portable, low-cost RNA sequencing of single cells at high of a drug targeting an intrinsically disordered protein involved in pancreatic throughput. Nat Methods. 2017;14:1–8. adenocarcinoma. Sci Rep. 2017;7:39732. 29. Yin X, Farin HF, van Es JH, Clevers H, Langer R, Karp JM. Niche-independent 52. Gulati AS, Shanahan MT, Arthur JC, Grossniklaus E, von Furstenberg RJ, high-purity cultures of Lgr5+ intestinal stem cells and their progeny. Nat Kreuk L, et al. Mouse background strain profoundly influences Paneth Methods. 2014;11:106–12. cell function and intestinal microbial composition. PLoS One. 2012;7: e32403. 30. Yin X, Mead BE, Safaee H, Langer R, Karp JM, Levy O. Engineering stem cell organoids. Cell Stem Cell. 2016;18:25–38. 53. Bevins CL, Salzman NH. Paneth cells, antimicrobial peptides and maintenance of intestinal homeostasis. Nat Rev Microbiol. 31. McLean WJ, Yin X, Lu L, Lenz DR, McLean D, Langer R, et al. Clonal expansion of Lgr5-positive cells from mammalian cochlea and high-purity 2011;9:356–68. generation of sensory hair cells. Cell Rep. 2017;18:1917–29. 54. Zhang Q, Pan Y, Yan R, Zeng B, Wang H, Zhang X, et al. Commensal 32. van Es JH, Jay P, Gregorieff A, van Gijn ME, Jonkheer S, Hatzis P, et al. Wnt bacteria direct selective cargo sorting to promote symbiosis. Nat Immunol. signalling induces maturation of Paneth cells in intestinal crypts. Nat Cell 2015;16(9):918–26. Biol. 2005;7:381–6. 55. Cunliffe RN, Rose FR, Keyte J, Abberley L, Chan WC, Mahida YR. Human 33. VanDussen KL, Carulli AJ, Keeley TM, Patel SR, Puthoff BJ, Magness ST, et al. defensin 5 is stored in precursor form in normal Paneth cells and is Notch signaling modulates proliferation and differentiation of intestinal expressed by some villous epithelial cells and by metaplastic Paneth cells in crypt base columnar stem cells. Development. 2012;139:488–97. the colon in inflammatory bowel disease. Gut. 2001;48:176–85. 56. Beumer J, Clevers H. Regulation and plasticity of intestinal stem cells during 34. Tian H, Biehs B, Chiu C, Siebel CW, Wu Y, Costa M, et al. Opposing activities homeostasis and regeneration. Development. 2016;143:3639–49. of notch and wnt signaling regulate intestinal stem cells and gut homeostasis. Cell Rep. 2015;11:33–42. 57. Janda CY, Dang LT, You C, Chang J, de Lau W, Zhong ZA, et al. Surrogate 35. Buczacki SJ, Zecchini HI, Nicholson AM, Russell R, Vermeulen L, Kemp R, Wnt agonists that phenocopy canonical Wnt and β-catenin signalling. et al. Intestinal label-retaining cells are secretory precursors expressing Lgr5. Nature. 2017;545(7653):234–7. Nature. 2013;495:65–9. 58. Gjorevski N, Sachs N, Manfrin A, Giger S, Bragina ME,Ordóñez-Morán P et al. 36. von Furstenberg RJ, Gulati AS, Baxi A, Doherty JM, Stappenbeck TS, Gracz Designer matrices for intestinal stem cell and organoid culture. Nature 2016; AD, et al. Sorting mouse jejunal epithelial cells with CD24 yields a 539:560–564. https://doi.org/10.1038/nature20168. population with characteristics of intestinal stem cells. AJP Gastrointest Liver 59. Sato T, Vries RG, Snippert HJ, van de Wetering M, Barker N, Stange DE, et al. Physiol. 2011;300:G409–17. Single Lgr5 stem cells build crypt-villus structures in vitro without a 37. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. mesenchymal niche. Nature. 2009;459:262–5. Gene set enrichment analysis: A knowledge-based approach for interpreting 60. Elias JE, Gygi SP. In: Hubbard SJ, Jones AR, editors. Target-Decoy Search genome-wide expression profiles. Proc Natl Acad Sci. 2005;102:15545–50. Strategy for Mass Spectrometry-Based Proteomics. Totoewa: Humana Press; 2010. p. 55–71. http://link.springer.com/10.1007/978-1-60761-444-9. 38. Mootha VK, Lindgren CM, Eriksson K-F, Subramanian A, Sihag S, Lehar J, et al. PGC-1α-responsive genes involved in oxidative phosphorylation are 61. Nesvizhskii AI, Aebersold R. Interpretation of Shotgun proteomic data. Mol coordinately downregulated in human diabetes. Nat Genet. 2003;34:267–73. Cell Proteomics. 2005;4:1419–40. 39. Xie X, Lu J, Kulbokas EJ, Golub TR, Mootha V, Lindblad-Toh K, et al. 62. Phanstiel DH, Brumbaugh J, Wenger CD, Tian S, Probasco MD, Bailey DJ, Systematic discovery of regulatory motifs in human promoters and 3’ UTRs et al. Proteomic and phosphoproteomic comparison of human ES and iPS by comparison of several mammals. Nature. 2005;434:338–45. cells. Nat Methods. 2011;8:821–7. Mead et al. BMC Biology (2018) 16:62 Page 24 of 24 63. Smyth GK. Linear Models and Empirical Bayes Methods for Assessing Differential Expression in Microarray Experiments. Stat Appl Genet Mol Biol. 2004;3:1–26. 64. Benajmini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc B. 1995;57:289–300. 65. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4:44–57. 66. Huang DW, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: Paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009;37:1–13. 67. Macosko EZ, Basu A, Satija R, Nemesh J, Shekhar K, Goldman M, et al. Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets. Cell. 2015;161:1202–14. 68. Satija R, Farrell JA, Gennert D, Schier AF, Regev A. Spatial reconstruction of single-cell gene expression data. Nat Biotechnol. 2015;33:495–502. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png BMC Biology Springer Journals
Loading next page...
 
/lp/springer_journal/harnessing-single-cell-genomics-to-improve-the-physiological-fidelity-q0XM7N7KLu
Publisher
BioMed Central
Copyright
Copyright © 2018 by Karp et al.
Subject
Life Sciences; Life Sciences, general
eISSN
1741-7007
D.O.I.
10.1186/s12915-018-0527-2
Publisher site
See Article on Publisher Site

Abstract

Background: Single-cell genomic methods now provide unprecedented resolution for characterizing the component cell types and states of tissues such as the epithelial subsets of the gastrointestinal tract. Nevertheless, functional studies of these subsets at scale require faithful in vitro models of identified in vivo biology. While intestinal organoids have been invaluable in providing mechanistic insights in vitro, the extent to which organoid-derived cell types recapitulate their in vivo counterparts remains formally untested, with no systematic approach for improving model fidelity. Results: Here, we present a generally applicable framework that utilizes massively parallel single-cell RNA-seq to compare cell types and states found in vivo to those of in vitro models such as organoids. Furthermore, we leverage identified discrepancies to improve model fidelity. Using the Paneth cell (PC), which supports the stem cell niche and produces the largest diversity of antimicrobials in the small intestine, as an exemplar, we uncover fundamental gene expression differences in lineage-defining genes between in vivo PCs and those of the current in vitro organoid model. With this information, we nominate a molecular intervention to rationally improve the physiological fidelity of our in vitro PCs. We then perform transcriptomic, cytometric, morphologic and proteomic characterization, and demonstrate functional (antimicrobial activity, niche support) improvements in PC physiology. Conclusions: Our systematic approach provides a simple workflow for identifying the limitations of in vitro models and enhancing their physiological fidelity. Using adult stem cell-derived PCs within intestinal organoids as a model system, we successfully benchmark organoid representation, relative to that in vivo, of a specialized cell type and use this comparison to generate a functionally improved in vitro PC population. We predict that the generation of rationally improved cellular models will facilitate mechanistic exploration of specific disease-associated genes in their respective cell types. Keywords: Single-cell RNA-seq, Chemical biology, Stem cell-derived models, Paneth cell, Intestinal organoid, Intestinal stem cell, Differentiation, Systems biology * Correspondence: bemead@mit.edu; jmkarp@bwh.harvard.edu Benjamin E. Mead and Jose Ordovas-Montanes contributed equally to the work as co-first authors. Alexandra P. Braun and Lauren E. Levy contributed equally to the work as co-second authors. Division of Engineering in Medicine, Department of Medicine, Brigham & Women’s Hospital, Harvard Medical School, Boston, MA, USA Full list of author information is available at the end of the article © Karp et al. 2018 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Mead et al. BMC Biology (2018) 16:62 Page 2 of 24 Background and have been identified in animal models to cause clear Intestinal organoids, derived from intestinal stem cells PC phenotypes, including lower DEFA expression [20], (ISCs) and composed of ISCs, Paneth cells (PCs), enter- defects in autophagy, granule formation, and secretion oendocrine cells (EECs), goblet cells, and absorptive [19, 21], and uncompensated ER stress [22]. While in vivo enterocytes, have been invaluable to the study of intes- models currently provide the most physiologically repre- tinal biology [1]. Recent advances in massively parallel sentative system to probe PC biology, they are inherently single-cell RNA-sequencing (scRNA-seq) [2] have en- complex and poorly scaled, hindering basic research into abled the cataloging of cell types and states of the molecular mechanisms of disease and the potential for murine small intestinal epithelium [3] and intestinal scalable therapeutic screening. Recently, conventional organoids [4], offering extensive insight into tissue het- intestinal organoids were used to describe the dynamics of erogeneity, specifically within subsets of rare secretory PC degranulation in response to multiple agonists [23] cell populations. However, there have been no formal and to assess PC suppression of enteric pathogens [24]. comparisons of how the in vitro intestinal organoid While these organoid studies are arguably more repre- condition recapitulates the defined in vivo cell types. sentative than other in vitro systems, the question of While the generation of comprehensive cellular atlases physiological fidelity of this heterogeneous system re- has become a major focus of a global effort to map mains unanswered, especially given that the timescales tissues in humans, model organisms, and derived to derive conventional organoids is typically less than a organoids at single-cell resolution [5], the challenge of week, while in vivo the lifespan of a PC is on the order how to functionally investigate key insights from cell of several weeks. types in vivo, or even more simply to confirm the To improve the representation of specific cell types in high-fidelity representation of these states in existing intestinal organoids, investigators have utilized cellular model systems remains [6, 7]. engineering approaches starting with ISCs to derive Intestinal organoids are a compelling system with multiple enriched or specialized models. These include which to study specialized cells of the epithelium. They enterocytes with improved intestinal ion transport [25], are self-organizing, stem cell-derived structures, which, epithelial monolayers capable of secretion and IgA to a reasonable degree, resemble their in vivo counter- transcytosis [26], and organoids enriched for the rare part, can be rapidly grown, and are amenable to many secretory EEC population [27]. However, in each in- biochemical and genetic perturbations [8]. Recent work stance, there has been no global comparison of the has demonstrated the utility of organoids in assessing extent to which intestinal organoids, or further special- bulk phenotypes that are readily observed and easily ized derivatives, recapitulate defined in vivo cell types selected for, such as phenotypes of cystic fibrosis and the and states. Moving beyond the generation of in vivo study of cancer-associated mutational signatures [9, 10]. tissue maps towards mechanistic insights, particularly in However, the application of organoid models to the study disease settings, will require an understanding of how of complex disease, such as polygenic inflammatory the in vitro organoid models utilized for such studies disease, has been limited. In such instances, the subtler represent the cell types and states identified beyond phenotypes, such as those present in inflammatory bowel single marker genes. disease (IBD), may not manifest if the originating cell state Here,weprovideaglobalcomparisonbetween thein present in vivo is not accurately represented within an vivo cellstatesofthe murine smallintestinalepithelium organoid model. This challenge is particularly clear in IBD and the in vitro conventional intestinal organoid, and [11], where loci identified through genome wide associ- establish a systematic workflow for improving the physio- ation studies (GWAS) have proven difficult to efficiently logical representation of stem cell-derived cell states to examine through the use of in vivo animal models. enable the creation of high-fidelity in vitro models. Taking For instance, PC dysfunction is implicated in Crohn’s the PC as a test case, we utilize massively parallel single-cell disease, a subset of IBD typically afflicting the small bowel transcriptomics (scRNA-seq; Seq-Well) [28] to bench- [12]. Co-localized with LGR5 ISCs of the small intestinal mark the conventional organoid model against its in vivo crypts, long-lived PCs [13] support maintenance of the counterpart and identify differences in developmental ISC niche, producing the Wnt and Notch signaling ligands pathway signaling between in vitro and in vivo cell states. WNT3 and DLL4 [14], and are potent modulators of the Single-cell transcriptomic approaches were key in enab- gut microflora through secretion of multiple antimicro- ling this study as epithelial cell types are challenging to re- bials including lysozyme (LYZ), phospholipase A2 group liably and prospectively isolate by fluorescence-activated 1B (PLA2G1B), angiogenin ribonuclease A family member cell sorting (FACS) due to the absence of robust surface 5 (ANG5), and alpha-defensins (DEFAs), amongst others markers and the spectrum of differentiation states present. [15]. Multiple allelic variants of NOD2, ATG16L1,and This profiling guides the rational augmentation of signal- XBP1, are associated with ileal Crohn’s disease [16–19] ing pathway activity during stem cell differentiation with a Mead et al. BMC Biology (2018) 16:62 Page 3 of 24 small molecule chemical induction method we previously our study to relate organoid-derived cell states to in vivo validated to enhance global Lyz gene expression [29]. We PCs. They are fully inclusive of the 14 high confidence validate our approach by generating an enhanced in vitro markers described for Paneth cells from the terminal physiological mimic of the in vivo PC and provide a ileum in the recently published mouse small intestinal detailed characterization of the derived cell state through atlas [3]. Of note, we extended our gene list beyond truly morphologic, proteomic, transcriptomic, and functional specific marker genes that are not expressed in other cell assays based on known signatures of in vivo PCs. Further- types as we were interested in a more comprehensive set more, we use our enhanced model and findings from its of PC-enriched genes for further comparison. transcriptomic and proteomic characterization to identify We next performed scRNA-seq using Seq-Well on Nupr1 as a potential stress-response factor that facilitates conventional organoids derived from a single donor the survival of PCs, demonstrating the improved ability ISC-enriched state (Fig. 1a). Beginning with murine small to examine gene function in vitro within a more repre- intestinal crypts, we directly enriched for LGR5+ ISCs sentative cell type. over 6 days following isolation within a Matrigel scaffold and medium containing recombinant growth factors EGF Results (E), Noggin (N), and R-spondin 1 (R), small molecules Using the PC to benchmark cell type representation of CHIR99021 (C), and valproic acid (V), as well as Y-27632 conventional organoids against their in vivo counterparts for the first 2 days to inhibit rho kinase and mitigate Conventional intestinal organoids produced from the anoikis, as previously described (ENR+CV) [29]. To spontaneous differentiation of ISCs have been used to ensure reproducibility within our system and limit the risk study PCs in vitro in multiple contexts [23, 24]. These in of interference in our chemical induction approach, we vitro PCs exist as part of a heterogeneous system, yet to conducted our study exclusively with recombinant growth be rigorously benchmarked against their in vivo counter- factors and not cell line-derived conditioned media. Cells parts. To better understand the composition of PCs were passaged into conventional ENR culture for an within conventional organoids and how well those PCs additional 6 days to allow multi-lineage differentiation approximate their in vivo counterparts, we sought to and produce stem cell-derived in vitro PCs. Following globally compare the conventional organoid-derived PCs scRNA-seq, we computationally identified six clusters and their in vivo counterparts through a single-cell tran- (amongst 2513 cells × 16,198 genes meeting quality scriptomic approach (Fig. 1a). standards, see Methods) in ENR organoids, which we To relate the organoid-derived PC state to in vivo PCs, label as ENR1-4, and EEC-1 and -2 for two EEC we first generated an unbiased reference in vivo scRNA-seq types (Fig. 1d). We identified ENR-4 as the cluster data set. We performed massively parallel scRNA-seq using most enriched for Lyz1 and our PC reference gene set the recently developed Seq-Well platform [28] on epithelial (effect size 0.721, ENR-4 vs. all ENR, *t test p <2.2 × −16 cells from the ileal region of the small intestine acquired as 10 ; for effect size details see Methods)(Fig. 1e, f). two biological replicates (see Methods). We assessed quality Having identified ENR-4 as the cell state of interest metrics for the number of genes, unique molecular identi- in organoids, we directly compared the top 200 most fiers (UMIs), mitochondrial genes, and ribosomal genes, all PC-like cells in ENR-4 to in vivo PCs by performing of which fell within expectations (all cells average: 1043 differential expression analysis (Fig. 1g). In comparing genes, 2168 UMIs, 5.4% ribosomal genes, 10.4% mito- the two cell types, it became evident that the majority of chondrial genes). UMI-collapsed cell-by-gene (7667 genes enriched by in vivo PCs were defensins and antimi- cells × 17,505 genes) expression matrices were analyzed crobials, including Defa22, Defa21, Zg16, Ang4, Defa3, −37 using Seurat (see Methods), performing dimensionality and Lyz1 (all p <2.92×10 , bimodal test, Bonferroni reduction, graph-based clustering, and deriving lists of corrected for multiple comparisons) (Fig. 1g, h). ENR-4 cluster-specific genes in order to identify PCs. Within cells were enriched for Chgb, an enteroendocrine the spectrum of cell types, we identified two clusters (2 marker, and translational biosynthetic genes likely and 11) enriched for Lyz1 expression (Fig. 1b, c), of indicative of the high rates of proliferation present which we determined cluster 11 to be fully mature PCs in ENR organoids (Fig. 1g). We further note the (n = 189 cells) based on uniform expression of a set of difference in genes arising from non-sex matched associated antimicrobial peptide marker genes such as comparison, like Xist, as a limitation of our compari- Defa22, Defa21,and Ang4 (receiver operating charac- son between a single donor for organoid derivation. teristic (ROC) test, area under the curve (AUC) > 0.99 Beyond these selected genes, we note a global reduc- for markers listed; cluster 11 average: 866 genes, 3357 tion in the fraction of the transcriptome of ENR-4 UMI, 3.5% ribosomal genes, 4.8% mitochondrial genes) cells producing the total cadre of in vivo PC marker (Additional file 1: Table S1). We further utilized these genes (effect size 1.25, InVivo vs. ENR, *t test p < −16 genes (genes with AUC > 0.65 for in vivo PC) throughout 2.2 × 10 ), suggesting that the current in vitro Mead et al. BMC Biology (2018) 16:62 Page 4 of 24 ab c d ef gh Fig. 1 (See legend on next page.) Mead et al. BMC Biology (2018) 16:62 Page 5 of 24 (See figure on previous page.) Fig. 1 Transcriptional benchmarking of in vitro Paneth cells (PCs) to in vivo. a Schematic of intestinal epithelial cell isolation from terminal ileum for unbiased identification of in vivo PC signature genes, and system for intestinal stem cell (ISC) enrichment to characterize in vitro PCs, via high-throughput scRNA-seq. b Marker gene overlay for binned count-based expression level (log(scaled UMI + 1)) of Lyz1, a canonical PC marker gene, on a tSNE (t-stochastic neighbor embedding) plot of 7667 small intestinal epithelial cells isolated from the terminal ileum; receiver operating characteristic (ROC)-test area under the curve (AUC) = 0.995, n = 2 mice, independent experiments (Additional file 1:Table S1). c Violin plot for the count-based expression level (log(scaled UMI + 1)) of Lyz1 across clusters identified through shared nearest neighbor (SNN) analysis (see Methods) over small intestinal epithelial cells; n = 196 cells in cluster 11, 7667 cells in total. d A tSNE plot of 2513 cells, with clusters identified through SNN (Additional file 1: Table S1 for full gene lists with ROC-test AUC > 0.60) from conventional ENR organoids; n =6 wells of ENR organoids. e Marker gene overlay for binned count-based expression level (log(scaled UMI + 1)) of Lyz1 on a tSNE plot from; ROC-test AUC = 0.856. f Violin plot of expression contribution to a cell’s transcriptome of PC genes across ENR organoid clusters from (d) (In vivo PC gene list −16 AUC > 0.65, Additional file 1: Table S1); effect size 0.721, ENR-4 vs. all ENR, *t test p <2.2 ×10 . g Row-normalized heatmap of top differentially expressed genes using bimodal test over single-cells from the top 200 PC-like cells from ENR-4 and the 196 in vivo PCs (cluster 11, from (c)); *bimodal −16 test, all displayed genes p <1.89 ×10 or less with Bonferroni correction. h Violin plots for the count-based expression level (log(scaled UMI + 1)) of −37 Lyz1, Ang4,and Defa3 in ENR and in vivo PCs; *bimodal test, all p <2.92 ×10 or less with Bonferroni correction. i Violin plot of expression contribution −16 to a cell’s transcriptome of PC genes (effect size 1.25, InVivo vs. ENR, *t test p <2.2 × 10 ), Wnt pathway (effect size 0.559, InVivo vs. ENR, *t test −8 −7 p <2.035 ×10 ) and Notch pathway (effect size −0.500, InVivo vs. ENR, *t test p <5.25 ×10 ) genes (see Additional file 2: Table S2 for gene lists) organoid-derived PCs are suboptimal for physio- plateaued by day 6 in ENR+CD versus ENR popula- logical studies (Fig. 1i). tions. Lgr5 expression in ENR+CD at 2 days versus Modulating key developmental pathways of stem ENR showed an insignificant plateau of expression, cell-derived systems has emerged as a paradigm in which trended down by 6 days. This may be indicative bioengineering to rationally generate cell types for basic of an expansion in ‘label-retaining’ secretory precur- research and therapeutic aims [30, 31]. Specifically, modu- sors [35]. The precursor population ENR + CV had no lating Wnt and Notch signaling has been suggested in the significant difference in PC or ISC markers relative to literature to increase the frequency and magnitude of Lyz1 ENR. ThesignificantincreaseinPCgeneexpression expression and protein in ISC-derived cells [29, 32–34]. in ENR + CD relative to ENR and ENR+CV over the Leveraging the single-cell transcriptomes of our in 6-day treatment suggests rapid enrichment following vitro- and in vivo-derived PCs, we confirmed that CI, supporting our hypothesis that alterations in Wnt Wnt target genes are enriched in vivo relative to in and Notch result in superior PC enrichment in vitro. vitro PCs (effect size 0.559, InVivo vs. ENR, *t test To phenotypically describe PC enrichment following −8 p <2.035 ×10 ) and Notch target genes were decreased CI, we performed imaging and immunocytochemistry −7 (effect size −0.500, InVivo vs. ENR, *t test p <5.25×10 ) for PC-associated features. After 6 days of ENR + CD, (Fig. 1i, Additional file 2: Table S2). As a result, we sought cell populations exhibited darkened annular morphology to comprehensively test if driving Wnt and inhibiting consistent with increased numbers of granule-rich cells Notch truly results in a more physiologically representa- (Additional file 3: Figure S1A). Confocal microscopy of tive PC versus the organoid-derived PC, beyond bulk mea- whole cell clusters stained for anti-DEFA and anti-LYZ sures of increased Lyz1 expression. showed an increase in LYZ+ and DEFA+ cells in ENR + CD compared to both ENR and ENR + CV (Fig. 2c). Rationally guided chemical induction of Wnt and Single-cell counting of confocal imaging showed a inhibition of Notch drives PC marker enrichment significant increase of DEFA and LYZ co-staining cells Beginning with an LGR5+ ISC-enriched population (ENR in ENR + CD (20–30% of cells) versus either ENR or +CV), we sought to profile how the modulation of Wnt ENR + CV (both < 5%; adj. p = 0.0001) (Additional file 3: and Notch signaling through small molecule inhibitors Figure S1B). Additionally, normalized z-axis profiles of in- would alter the in vitro PC state, as suggested by our dividual co-staining cells within cell clusters revealed a con- transcriptional profiling. We performed chemical induction sistent distribution of DEFA (luminally polarized) and LYZ (CI) using the previously identified compounds C to drive (diffuse) (Additional file 3:FigureS1C 1–3). High-resolution Wnt signaling and DAPT (D), a gamma-secretase inhibitor, fluorescent imaging of individual co-staining cells from to inhibit Notch (ENR+CD) (Fig. 2a)and measured bulk freshly isolated small intestinal crypts (in vivo equivalent) gene expression of PC (Lyz1, Defa1, Mmp7)and ISC(Lgr5) and 6-day ENR + CD-treated cells showed a similar markers every 2 days for a total of 6 days (Fig. 2b). ENR polarized distribution of LYZ- and DEFA-stained gran- +CD-treated cells had statistically significant increases ules, although freshly isolated cells appeared to be in Lyz1 (adj. p = 0.005, see Methods)and Mmp7 (adj. more granular than CI-PCs (Fig. 2d). p = 0.005) within 2 days compared to ENR, with differ- To confirm the extent of enrichment seen in whole ences plateauing around 4 days. Defa1 (adj. p = 0.004) population imaging, the prevalence of PCs in ENR + CD expression was significantly increased by day 4 and relative to ENR was assessed by flow cytometry over the Mead et al. BMC Biology (2018) 16:62 Page 6 of 24 ab c de f gh Fig. 2 (See legend on next page.) Mead et al. BMC Biology (2018) 16:62 Page 7 of 24 (See figure on previous page.) Fig. 2 Establishing chemically induced Paneth cell (PC)-enriched cultures. a Schematic of small molecule-driven differentiation of LGR5+ ISCs (C - CHIR99021, D - DAPT) and non-specific differentiation. b mRNA expression of PC (Lyz1, Defa1, Mmp7) and ISC (Lgr5) markers relative to ENR, for ENR+CV and ENR + CD at 2 (D2), 4 (D4), and 6 days (D6) (n = 3 biological replicates; two-way ANOVA with multiple comparison test vs. ENR; ** adj. p < 0.01, *** adj. p < 0.001). c Representative confocal imaging of whole cell clusters for PC antimicrobials following 6 days in ENR+CD versus ENR and ENR+CV: stained for anti-DEFA, anti-LYZ and counterstained with DAPI and for actin (phalloidin). d High-resolution fluorescent imaging of in vivo and in vitro single cells from 6-day culture in ENR + CD shows similar morphology and antimicrobial expression: stained for DEFA and LYZ, and counterstained with DAPI and for actin (phalloidin). e Viable cell populations from ENR, ENR+CD, and ENR+CV precursor culture have distinct populations based on CD24 and LYZ content, indicative of PC maturity (n = 3 biological replicates; ENR + CV, days 4, 6, 12, n = 2 biological replicates day 8). f Volcano plot of differentially regulated proteins between 6-day (6D) ENR + CD and ENR cells shows clear enrichment in secreted and PC-associated proteins (labeled). Cut-offs are 2 standard deviations outside the mean expression level of the set and FDR < 0.05. g Rank-order log fold change of detected PC antimicrobial proteins and between 6-day ENR + CD and ENR cultures (n =4). h Rank-order log fold change of detected secretory proteins associated with EEC and goblet lineages in ENR + CD relative to ENR cultures (n =4) course of 12 days, a longer term culture than typical for samples and analyzed them in a single 10-plex by conventional organoids. We identified an in vivo PC LC-MS/MS (Additional file 4: Figure S2A). We identi- phenotype as CD24 and LYZ co-positive cells, as per fied 8015 unique proteins within all samples; each previous reports [36], and noted the presence of replicate pair (ENR + CD/ENR) was normally distributed single-positive LYZ+ or single-positive CD24+ popula- (not shown) and correlated with all others, indicating con- tions, indicative of alternative cell differentiation, imma- sistentproteomeenrichment(Additional file 4:FigureS2B). ture, or non-physiological PCs (representative populations We looked at the sample pairs in aggregate and classified Additional file 3:FigureS1D,representative gating proteins significantly enriched in ENR + CD and ENR by a Additional file 3: Figure S1E). ENR + CD had substantial false discovery rate (FDR) < 0.05 and log fold change (± 2σ) enrichment at all time points for double-positive, and (Fig. 2f and Additional file 5: Table S3). There were 249 single-positive LYZ+ or CD24+ populations relative to ENR + CD-enriched proteins, 212 ENR-enriched proteins, ENR, as well as a consistent decrease in the double and 7553 shared proteins. Known PC markers, including negative population in agreement with the PC phenotype LYZ, DEFAs, and other secretory pathway components, (Fig. 2e). Notably, both ENR and ENR + CD experience were identified as significantly enriched in ENR + CD ver- declines in total cell viability, with ENR + CD having sus ENR alone. Of known antimicrobial proteins produced greater survival at longer times, suggesting both a reduc- by PCs, we detected 10 DEFAs, 5 CRS peptides, 6 ribonu- tion in anoikis, a potentially physiological ‘long-lived’ PC cleases, 12 lectins, LYZ1, and PLA2G1B with differential phenotype in ENR + CD versus ENR, or an enhance- abundance between ENR + CD and ENR (Fig. 2g). Each ment in niche-supporting functionality (Additional class of antimicrobials had at least one ENR + CD enriched file 3: Figure S1F). Overall, imaging and flow cytometry protein (+ 2σ), with the ribonucleases significantly enriched demonstrate a significant increase in cells morphologically and a majority of the lectins and DEFAs unregulated be- resembling in vivo PCs with respect to granularity, polar- tween the two conditions. Proteins associated with the EEC ity, and antimicrobial co-expression in ENR + CD com- lineage (secretogranins, chromogranins, and neuropeptides) pared to conventional ENR organoids (Fig. 2c–e and were also enriched in ENR + CD, in addition to multiple Additional file 3:FigureS1A–F). other secreted components, including Wnt ligands, and the complement pathway components C3 and CFI (Fig. 2h). In Chemically induced PC proteome is enriched for sum, we see a broad diversity of PC-associated antimicro- components of secretory lineages bials with some enrichment of EEC-associated proteins in With ENR + CD apparently providing a more prevalent ENR + CD relative to ENR. and physiological PC population, we sought to more Additionally, we characterized enriched biological func- globally characterize the differences between in vitro tions, cellular compartments, and molecular functions PCs (ENR vs. ENR + CD) at 6 days. Because our single using DAVID v6.8 and the gene ontology database. All cell transcriptomic comparison revealed that many of sets had high database coverage (greater than 85%) of the differential genes between PCs in conventional orga- queried proteins. The ENR + CD proteome is signifi- noids and in vivo were lineage-defining protein prod- cantly enriched for extracellular and protein process- ucts, we sought to assess the total intracellular proteome ing compartments and secretory-associated functions between the conventional organoid and our chemically (Additional file 4: Figure S2C), while the ENR prote- induced model through liquid chromatography mass ome favors translation, intracellular compartments, spectrometry (LC-MS/MS)-based proteomics. We quan- and translational activities (Additional file 4:Figure tified relative protein abundance using isobaric mass S2D). Of note, there are the extracellular exosome and tag labeling from four ENR and four ENR + CD calcium ion-binding associated proteins in the ENR + Mead et al. BMC Biology (2018) 16:62 Page 8 of 24 CD proteome that are indicative of the intestinal cells expressing antimicrobial Lyz1, Defa24, Defa3, Mmp7, epithelial secretory phenotype (for a complete list of and EEC marker Chga, and that ENR enriched for absorp- DAVID enrichments, refer to Additional file 6:Table tive marker Fabp2-expressing cells (Fig. 3b). S4). These functional enrichments further support the To assess subpopulation structure and provide a notion that the ENR + CD-cultured organoids are more robust measure of composition beyond canon- enriched in secretory cells, including PCs, although it ical marker genes, we performed unsupervised KNN does not rule out potential co-enrichment for the EEC graph-based clustering on the captured cells (Fig. 3c, d lineage. Finally, we sought to identify transcription and Additional file 1: Table S1 for full gene lists), factors (TFs) that may mediate PC-specific differenti- distinguishing four clusters in each treatment condi- ation using GSEA [37, 38]withthe MSigDB transcrip- tion. We then scored individual clusters according to tion factor target (v5.2) gene set database [39]witha the amount of the transcriptome within each cell dedi- moderately conservative cutoff (see Methods). We cated to synthesizing the respective enriched proteins generated an enrichment map [40, 41]ofseveral TF from the bulk proteome data. We observed that ENR targets significantly enriched in both the ENR + CD + CD clusters yield a significant enrichment for those and ENR proteomes. In ENR + CD, the nuclear recep- proteins detected in the up-regulated proteome (effect −16 tors for progesterone, aldosterone, and glucocorticoid, size 1.38 ENR + CD vs. ENR clusters, p <2.2 ×10 ) as well as the cellular differentiation-implicated TAL1, and that the down-regulated proteins were enriched in RP58, and NRSF, were significantly enriched. In ENR, the ENR and ENR + CV conditions (Fig. 3d, e and data the primary known enrichment was for the cell cycle and not shown). Intriguingly, at the level of clusters, the proliferation-related E2F TF family (Additional file 4: upregulated proteome was not evenly distributed Figure S2E). These potential TFs are consistent with across all cells in ENR + CD, but was most enriched in CI-PC treatment driving expected terminal differentiation cluster ENR + CD-4, which represented approximately of specialized cells, as opposed to conventional organoid 10% of ENR + CD cells (effect size 2.40 ENR + CD-4 −16 culture, which supports a broad mix of intestinal epithelial vs. all cells, p <2.2 ×10 )(Fig. 3d, e). cells, including proliferating populations. Furthermore, To address ENR + CD composition and how it relates to this analysis suggests potential targets, such as progester- conventional organoids, we interrogated the expression of one, aldosterone, and glucocorticoid, to modulate the Lyz1, Chga, and other selected genes across each cluster differentiation programs of this secretory cell population (Fig. 4a). We noted that clusters ENR-4 and ENR + CD-4 in future studies. shared expression of Lyz1, Defa24, Defa3,and Mmp7, yet ENR + CD-4 cells produced significantly more of each −74 Single-cell RNA sequencing profiles heterogeneity of canonical PC gene (bimodal test, p <6.80×10 for genes chemically induced PCs, revealing subsets with improved listed, Bonferroni corrected for multiple comparisons). transcriptional similarity to in vivo Furthermore, both ENR-4 and ENR + CD-4 cells lacked With the apparent co-enrichment of canonical PC and expression of EEC genes like Chga, which was observed in EEC proteins in the ENR + CD proteome, we sought to the EEC-1 and EEC-2 clusters arising from mixed-grouping identify whetherweproduceahomogenouspopulationof of the sample, as well as in ENR + CD-2 and ENR + CD-3 mixed-lineage secretory cells or a spectrum of unique cell (Fig. 4a). Altogether, this suggests that ENR + CD drives states between EEC and PC. We performed scRNA-seq PC differentiation while also inducing a secretory tran- using the Seq-Well platform on cells from ENR + CD and sition state (ENR + CD-2 and 3) expressing a mix of PC the precursor ENR + CV conditions to analyze alongside and EEC marker genes (Additional file 1:Table S1 for conventional ENR organoids. To ensure experimental full gene lists). robustness, we assessed quality metrics for the number of We next sought to compare the states generated in genes, UMIs, mitochondrial genes, and ribosomal genes by vitro to those observed in vivo with our refined system. cluster, all of which fell within expectation (Additional file 7: Using the gene list of in vivo PC markers and further de- Figure S3). UMI-collapsed digital gene expression matrices fining a list for in vivo EECs (see Methods) captured on were analyzed using Seurat (see Methods), and displaying the Seq-Well platform (Additional file 1: Table S1), we all three treatments (ENR + CV, ENR, ENR + CD) in tSNE observed that the percentage of a cell’stranscriptome space demonstrated clear separation between each condi- dedicated to synthesizing defining Paneth genes was tion (Fig. 3a), illustrating that the unique transcriptional significantly enriched relative to ENR-4 in clusters ENR −5 differences induced by each treatment were conserved + CD-2, -3, and -4 (effect size 0.15, p <3.43 × 10 ;ef- −16 across all cells. Plotting key genes demonstrated that, as fect size 0.829, p <2.2 ×10 ; effect size 2.52, p <2.2 × −16 expected, all cells expressed high levels of Epcam,that 10 , respectively) with an increase in expression of ENR + CV cells had enhanced Mki67, a marker of prolifer- EEC genes across ENR + CD-1, -2, and -3 but not ENR –16 ation, that the ENR + CD condition led to enrichment of + CD-4 (effect size 1.30, p <2.2 ×10 ; effect size 1.82, Mead et al. BMC Biology (2018) 16:62 Page 9 of 24 a b cd e Fig. 3 Single-cell RNA-sequencing reveals cellular composition across treatments and origins of proteomic data. a A tSNE plot of single cells derived from ENR + CV (n =985 cells), ENR (n =2544 cells), and ENR +CD (n = 2382 cells) harvested at day 6 of differentiation, colored by treatment; n =6 wells for each condition. b Marker gene overlays (on plot from (a)) for binned count-based expression level (log(scaled UMI + 1)) of individual genes of interest. c A tSNE plot, with clusters identified through SNN graph-based clustering (see Additional file 1: Table S1 for marker gene lists), highlighting distinct cell states within each organoid; opacity of density clouds correspond to the Paneth cell score of ENR-4, ENR + CD-3, and ENR + CD-4 clusters (see Fig. 4b). d Violin plot of expression contribution to a cell’s transcriptome of ENR + CD proteome-enriched genes across organoid clusters from (c) (see Additional −16 file 1:Table S1 forfullgene list);effect size 2.40 ENR+CD-4 vs.all cells, p <2.2 ×10 . e Frequency of each cluster observed within each organoid condition as a fraction of the total cells in each condition –16 −16 p <2.2 ×10 ; effect size 1.118, p <2.2 ×10 ; effect size ENR + CD-4 cells relative to in vivo PCs demonstrated a 0.0465, p = 0.2339, respectively) (Fig. 4b). Notably, ENR + striking similarity relative to the difference observed be- CD-4 cells (~10%) had a three-fold increase in the tran- tween in vivo and ENR-4 cells (PC fraction of in vivo scriptional resemblance to in vivo PCs relative to ENR-4 transcriptome: effect size 0.237 InVivo vs. ENR + CD-4, (53.4% of transcriptome ENR + CD-4 vs. 16.5% of tran- p < 0.0055; effect size 1.25 InVivo vs. ENR-4, p <2.2 × −16 scriptome ENR-4) (quantification of Fig. 4b). Furthermore, 10 (Additional file 1:Table S1). 45% of ENR + CD cells express a secretory PC-like tran- In Fig. 4c, we present a heatmap of scaled expression scriptional phenotype that is at least two-fold enhanced values for the top genes (AUC > 0.65) used for the in relative to conventional organoids (33.9% of transcriptome vivo Paneth score across ENR-4, ENR + CD-4, and the in ENR + CD-3 and -4 vs. 16.5% ENR-4). Comparing the vivo cluster used to define PCs. We observed that the Mead et al. BMC Biology (2018) 16:62 Page 10 of 24 Fig. 4 (See legend on next page.) Mead et al. BMC Biology (2018) 16:62 Page 11 of 24 (See figure on previous page.) Fig. 4 Transcriptional identity of chemically induced Paneth cells (CI-PCs) within conditions and related to in vivo PCs. a Violin plots for the count-based −74 expression level (log(scaled UMI + 1)) of selected genes across called clusters, colors correspond to clusters in Fig. 4c;*t test, p <6.80× 10 or less with Bonferroni correction, for Lyz1, Defa24, Defa3,and Mmp7 ENR + CD-4 relative to ENR-4. b Violin plot of expression contribution to a cell’s transcriptome of in vivo PC and enteroendocrine marker-cell genes (see Additional file 1: Table S1 for full gene list, AUC > 0.65); effect size 2.52 ENR + CD-4 vs. −16 ENR-4, p <2.2 × 10 for PC score; effect size 0.0465, p = 0.2339 ENR + CD-4 vs. ENR-4 for enteroendocrine cell score. c Row-clustered heatmap of z-scores (−2.5 to 2.5; purple to yellow) for defining genes (n = 69 with AUC > 0.65 of in vivo PCs, see Additional file 1:Table S1 forfullgene list) across top 200 cells for PC score (Fig. 5b) from ENR-4 and ENR + CD-4 conditions compared to two biological replicates of in vivo PCs from the terminal ileum (n = 196 cells) enhanced PC phenotype in ENR + CD-4 (effect size around 6 hours post-wash (two-way ANOVA, stimulant −16 1.144 ENR + CD-4 vs. ENR-4, p < 2.2 × 10 ) correlated p < 0.0001, time-point p < 0.0001) (Fig. 5b). The ob- with a greater expression of signature genes, such as served PC secretion in response to CCh is consistent Lyz1, Lyz2,and Defa5, and a greater diversity of anti- with observations made in ex vivo crypts, though over microbial peptide genes, such as Ang4, Defa3, and the appreciably longer time scales, likely due to the added metalloprotease Mmp7. diffusion barrier of the organoid structure and matrigel To confirm and extend our findings of pathway-based [44]. We next identified how LYZ secretion changes modulation, we scored clusters for enrichment or depletion over the course of differentiation. Beginning with an of canonical growth factor-induced pathways. CHIR99021 ISC-enriched population, we assayed for secreted LYZ activates the Wnt pathway, and we observed a significant in cell culture supernatants every 2 days for 6 days of enrichment for Wnt target genes in all CI-PC clusters ENR + CD culture, following a 24-h stimulation with −16 (effect size > 0.999, p < 2.2 × 10 for all ENR + CD clus- CCh or without (basal collection/non-stimulated). ters vs. ENR-4) (Additional file 8: Figure S4A). While Notable increases in functional secretion (stimulated DAPT is a Notch pathway inhibitor, levels of Notch target relative to basal) occurred at days 4 and 6 (two-way genes were largely greater than or equivalent to ENR-4 ANOVA, stimulant p < 0.0001, time-point p < 0.0001) cells across CI-PC clusters, except for significant depletion (Fig. 5a). Compared to conventional organoids and −16 in ENR + CD-4 (effect size −0.658, p <2.2 ×10 ENR + ISC-enriched precursors, ENR + CD secreted significantly CD-4 vs. ENR-4) (Additional file 8: Figure S4B). This more basal LYZ (p < 0.0001) and was the only population suggests that complete Notch suppression is key for PC that showed grossly measurable CCh-induced secretion differentiation distinct from an EEC fate. Additionally, (adj. p = 0.03) (Fig. 5c). This result is consistent with the given the recognized role for distinct respiratory potential observed enrichment between chemically induced popula- in enterocytes, ISCs, and PCs, we scored cells across tions relative to conventional. respiratory electron transport genes [42, 43]. ENR + CD-4 Based on the broad spectrum of antimicrobials de- had the lowest cluster score relative to all cell subsets tected proteomically, transcriptionally, and functionally, –16 (effect size −1.4649, p <2.2 ×10 )(Additional file 8: we hypothesized that ENR + CD possess greater bacteri- Figure S4C). Together, this suggests that Wnt signaling cidal effects than conventional organoids. We assayed is necessary but not sufficient to specify the mature PC for bacterial growth modulation by suspending cell clus- phenotype and that Notch and metabolic conditions ters with common laboratory strains of gram-negative play a larger role in the decision between PC and EEC and gram-positive bacteria in exponential growth. fates. CI-PCs significantly suppressed growth of gram-positive L. lactis MG1363 (adj. p = 0.0001), which did not occur Chemically induced PCs mimic in vivo stimulant-induced with conventional organoids, indicative of increased secretion and demonstrate selective modulation of PC-associated antimicrobial activity (Fig. 5d). Both ENR bacteria in co-culture (adj. p = 0.0005) and ENR + CD (adj. p = 0.01) co-culture In addition to our morphological, proteomic, and tran- showed significant increase in gram-negative E. coli scriptional characterization of PC phenotype in ENR + MG1655 growth but no appreciable effect was observed CD and ENR, we sought to measure physiological on the growth of gram-positive E. faecalis V583 versus function by assessing stimulant-induced secretion of bacteria alone (Fig. 5d). While this assay simplifies the antimicrobials. We assessed the dynamics of LYZ PCs’ physiological environment and may not be a direct accumulation in the supernatant media of cultures proxy for strain-specific growth modulation, it does following media wash, basally and after stimulation demonstrate that the PC-enrichment of ENR + CD ver- with carbachol (CCh), a cholinergic agonist known to sus conventional organoids enables detectable in vitro induce PC secretion [44]. CCh (10 μM) induced a rapid bacteria species-specific PC antimicrobial response, accumulation of LYZ within 2 hours that plateaued opening avenues for future experimentation. Mead et al. BMC Biology (2018) 16:62 Page 12 of 24 a bc de Fig. 5 Chemically induced Paneth cells (CI-PCs) are functional in response to host and microbial stimuli. a Supernatant LYZ from 24-h basal and 10 μm CCh-stimulated LYZ cells at varying number of days in ENR + CD culture (top). DNA content from matched samples (bottom) (n = 8 well replicates; SEM error bars too small to visualize). b Supernatant LYZ from 6-day ENR + CD collected basally and following 10 μm CCh stimulation for 0.5, 2, 4, 6, and 24 h (top). DNA content from matched samples basally and following 10 μm CCh stimulation (bottom) (n = 8 well replicates). c 24-h basal (non-stimulated) and 10 μm CCh-stimulated LYZ secretion in 6-day ENR + CD versus ENR and ENR + CV (n = 8 well replicates; two- way ANOVA with multiple comparison test; ns non-significant, * adj. p < 0.05, **** adj. p < 0.0001). d 4-h co-culture of freshly passaged 6-day ENR and ENR + CD cells and select gram-negative and gram-positive aerobic bacteria (n = 13 well replicates; two-way ANOVA with multiple comparison test, * adj. p < 0.05, *** adj. p < 0.001, **** adj. p < 0.0001). e Normalized cellular viability, caspase activity per viable cell, and cytotoxicity per viable cell from 24-h and 48-h ENR and ENR + CD co-cultures at specified mixing ratios (n = 9 well replicates from three biological donors; one sample t test,* p <0.05, ** p <0.01, *** p <0.001, **** p <0.0001) Chemically induced PCs provide niche support and co-culture viability, caspase activity, and cytotoxicity 24 enhance conventional organoid survival and 48 h following re-plating in ENR media. If there was Beyond the generation of antimicrobial peptides, PCs no appreciable interaction, positive or negative, between provide niche support for ISCs. We sought to test if the two populations we would expect to see a linear CI-PCs provided niche factors known to drive epithelial trend of every measured variable throughout mixing regenerative turnover. We performed co-culture experi- ratios. However, we observe a significant positive inter- ments, mixing and re-plating cell populations derived action where the presence of both populations drives an from 6 days of ENR or ENR + CD culture and assayed overall increase in cellular viability, beginning at 24 h Mead et al. BMC Biology (2018) 16:62 Page 13 of 24 (one-sample t test 1:1 p = 0.037) and increasing at 48 h for 2 days following a 6-day course of ENR + CD, where (one-sample t test 1:1 p =0.001 and 1:3 p < 0.001) again a profound, but not total, decline in cellular viability (Fig. 5e). This is likely due to a significant decrease in was observed. Further, it appears that TFP treatment is se- overall apoptosis relative to the total cell population lectively more toxic to PC and PC-progenitor populations (one-sample t test 24 h 1:1 p = 0.004 and 1:3 p = 0.032, relative to non-PC populations (Fig. 6d). In total, this ini- 48 h 1:3 p = 0.003), and unrelated to changes in cellular tial investigation suggests that NUPR1 may be a critical cytotoxicity. We believe that the presence of a TF in PC development and survival, which carries PC-enriched population (from ENR + CD) is driving therapeutic implications which we will seek to validate this effect by providing increased soluble regenerative with expanded gene-perturbation studies in vitro and in factors to the ISC population in ENR organoids, in- vivo in future work. creasing the generation of new cells, and resulting in a lower overall rate of apoptosis. Discussion We sought to directly compare a specific cell type present Mapping of in vivo PC-associated transcription factors to in vivo to that derived in an intestinal organoid in vitro, in vitro proteome and transcriptome reveals Nupr1 as with the main goal of understanding the nature and extent important in epithelial survival of divergence between the in vitro and in vivo conditions. Finally, we sought to use this physiologically improved in Empowered by recent advances in massively-parallel vitro PC system (ENR + CD) to identify novel factors po- scRNA-seq, we employed a generalizable approach to tentially supportive of PC survival or differentiation. Using define the relation between in vivo and organoid-derived our in vivo PC and EEC gene lists, and filtering for only in vitro PCs and employ a rationally identified interven- TFs (using TFdb, downloaded September 2017) [45], we tion to improve in vitro representation through chemical identified a set of PC- or EEC-specific TFs. We mapped modulation of developmental pathways. Our scRNA-seq these TFs to our in vitro proteome (Fig. 6a and Additional approach enabled the identification of populations of file 5: Table S3), which revealed the previously unreported interest both in vivo and in vitro and the analysis of differ- NUPR1 as the most enriched PC-specific TF in ENR + ences between subpopulations which would have other- CD. This finding was supported by differential expression wise been greatly obscured in a bulk analysis of this between ENR + CD-2 (most EEC-like cells) and ENR + heterogeneous system. While others have recently used −37 CD-4 (p <3.12×10 , bimodal test, Bonferroni corrected scRNA-seq to profile the heterogeneity of intestinal orga- for multiple comparisons) (Fig. 6b). We further identified noids [4] and the murine small intestinal epithelium [3], Nupr1 in our in vivo PC populations, which showed spe- we provide a direct comparison between the model and cific and enriched expression of Nupr1 by in vivo PCs system. This comparison is key to understanding what (ROC test, AUC = 0.833) (Fig. 6b). Intriguingly, Nupr1 is a complex models, such as intestinal organoids, really repre- stress-response gene, known to promote cellular survival sent, how they may best be utilized and rational strategies and senescence through mediation of autophagy, and has for improvement. primarily been studied in the context of cancer [46–48]. In our comparison, we identified that the PC-state of Autophagy and stress response have repeatedly been im- conventional intestinal organoids poorly represents the plicated through GWAS study in PCs in IBD; however, extent and breadth of antimicrobial gene expression, and Nupr1 has only ever been reported in a single IBD GWAS that modulation of Wnt and Notch during differentiation study, and its role in PC biology has not been formally in- may improve physiological representation. Using a com- vestigated [49]. With our model, we sought to test the role bination of small molecule promotion of Wnt and inhib- of NUPR1 on in vitro PC survival through the small ition of Notch signaling that had previously been shown to molecule inhibition of NUPR1 with trifluoperazine (TFP) improve the bulk expression of PC-marker gene Lyz1 [29], [50, 51]. While genetic perturbation may provide for more we drove a secretory differentiation program and enriched specific effect measurement, we chose to use TFP as a for mature PCs with greater diversity and expression of simple, albeit less specific, modulator, as the complexity antimicrobial peptides relative to existing in vitro models involved in selecting for a genetically perturbed popula- and, thus, are more representative of in vivo PCs. Imaging tion of PCs in organoids, if Nupr1 is a survival gene, is be- of this population revealed that they are positive for the yond the scope of the present study. We first tested how antimicrobials LYZ and DEFA, clearly polarized, and different dosages impact PC differentiation in combin- granulerich, suggestive of amaturePC. This population is ation with ENR + CD for 6 days, where doses above 1 μM approximately six-fold more abundant in ENR + CD than led to near total cell death, and where the few surviving an ENR organoid, as confirmed through image quantifica- cells were primarily non-PC (Fig. 6c). This suggests that tion, flow cytometry and scRNA-seq. We further character- Nupr1 is likely critical to cellular survival during the CI ized the subpopulation enrichments of our ENR + CD differentiation process. We also tested the addition of TFP culture and directly compared it to conventional organoids. Mead et al. BMC Biology (2018) 16:62 Page 14 of 24 ab cd Fig. 6 CI-PCs reveal putative function of Nupr1 transcription factor in Paneth cell (PC) survival. a ENR + CD is enriched for in vivo PC and EEC transcription factors, including Nupr1 (n = 4). b Violin plots for the count-based expression level (log(scaled UMI + 1)) of Nupr1 across in vivo and in vitro called clusters. c Nupr1 inhibition with trifluorperazine (TFP) treatment concurrent with 6-day ENR + CD differentiation reveals dose- dependent toxicity, with preference to PCs (CD24+ and LYZ+) and PC-like (CD24+ and LYZ+) populations as assessed by flow cytometry (n =3 biological replicates). d 2-day TFP treatment following 6-day ENR + CD differentiation reveals dose-dependent toxicity, with preference to PCs (CD24+ and LYZ+) and PC-like (CD24+ and LYZ+) populations as assessed by flow cytometry (n = 5 biological replicates) We identified two subpopulations in scRNA-seq (ENR + identify as the approximately 5% of single-staining LYZ+ CD-3 and ENR + CD-4) that account for approximately cells present in ENR organoids as assessed by flow cytome- half of the ENR + CD-treated cells with a high-degree of try (Fig. 2e). This finding makes sense, given that conven- transcriptional similarity to in vivo PCs, a greater percent- tional organoid culture often occurs on the time scale of a age/matching than the ENR-subpopulation that most re- week, while in vivo PCs are relatively long lived (several sembles an in vivo PC (ENR-4). From this analysis, we weeks), develop in non-sterile conditions, and presumably believe that in vitro PCs characterized in the past [23, 24] would require longer periods of culture to reach maturity likely represent secretory precursor populations lacking the in vitro. Indeed, our ENR + CD cultures show increasing full phenotypic repertoire of the in vivo PC, which we PC populations up to 12 days of culture and may likely Mead et al. BMC Biology (2018) 16:62 Page 15 of 24 continue gaining at longer time points. While our ap- inhibition is important for mature PC development, pos- proach moves us much closer to generating the in vivo PC sibly as a balance between differentiation and cell survival. (fraction of in vivo transcriptome: effect size 0.237 InVivo Interestingly, we also see a notable gradient in cellular vs. ENR + CD-4, p < 0.0055; effect size 1.25 InVivo vs. respiration across subpopulations, lowest in the PC and −16 ENR, p <2.2 ×10 ), we still do not capture the total highest in the stem cell and EEC lineages, in agreement amount of antimicrobial peptides present in vivo, and with recent work on the metabolic differences within the propose pathways to modulate in future studies. stem cell niche [43], as another potential cue to further Evidence suggests that PC antimicrobial expression specify PC differentiation. Future studies should incorpor- and function are influenced by genetic background [52] ate temporal aspects to growth factor delivery akin to and implicated in intestinal disease, including IBD [53]. what has been shown for degradable matrices [58]to How genetic background may influence differentiation enhance purity and yield, and explore the role of meta- through this protocol is yet to be studied but especially bolic utilization in addition to growth factor signaling in prudent, as we demonstrated the ability to detect a cellular fate determination. Indeed, while we implemented broad spectrum of antimicrobial proteins and peptides an approach of chemical induction to drive model im- and their differential abundance within a PC-enriched provement, there are many other approaches which may population. Interestingly, we identified that the same be implemented to seek similar shifts in the composition subpopulation (ENR + CD-4) with the most transcrip- of organoid systems, such as using novel ligands or tuning tional overlap to our bulk ENR + CD-enriched proteome of the material supports. Further, we appreciate that many also most closely resembles the in vivo PC. While this investigators have begun using R-spondin or Wnt condi- subpopulation does not account for the majority of ENR tioned medias in their intestinal organoid cultures; how- + CD-cultured cells, it appears that ENR + CD-4 consist- ever, to ensure consistency and control of the system ently drives the PC phenotype in vitro. In addition to upon chemical induction we chose to use recombinant assessing the role of genetic background or disease state growth factors in lieu of a less-characterized media prod- on antimicrobial content, our platform also affords the uct, but this may not preclude their usage. Overall, our ability to interrogate how alterations in protein process- analyses of single cell heterogeneity show that our sys- ing and storage in PCs affects the proteome, which has tem is well positioned to further investigate the effects been shown to drive shifts in the microbiome and may of both known and unknown physiological cues on PC be implicated in disease [54, 55]. Finally, while we differentiation and function. demonstrate an enriched phenotypic spectrum of anti- One of the most important features we established with microbials and Wnt ligands, we also identified several our CI-PCs was the ability to measure PC functional neuropeptides and hormone products associated with enrichment through simple soluble assays. We demon- the EEC lineage within our system. Given that multiple strated sufficient functional enrichment in PCs such that studies have linked the differentiation of PCs and EECs enzymatic activity assays can detect stimulant-induced through a common progenitor population [56], it is secretion of antimicrobials as well as the promotion of reasonable to expect enrichment in one population the ISC niche. Moreover, microbe co-culture assays would also allow for some overlap with the other, as we with our enriched cells produce measurable and select- see in our scRNA-seq data. Future experiments will seek ive microbial growth modulation not observed using to leverage these transition states to more formally iden- conventional organoids. Co-culture strains were chosen tify genetic programs that underlie common or unique to demonstrate proof of concept of selective antimicro- developmental trajectories. bial action and assess functionality compared to con- To understand how our chemical induction led to ventional organoids. Given the results showing selective distinct secretory subpopulations within the CI-PCs, we modulation of bacterial growth, we believe that our sys- mapped Wnt, Notch, and metabolic gene sets onto each tem could serve as a tool to further probe host-microbe subpopulation. In our system, Notch-signature is highest interaction in vitro. Furthermore, it would allow for in- in the stem cells and EECs, lower in enterocytes, and vestigations of both microbial mechanisms that elicit lowest in PCs. Our system’s Wnt signature is relatively PC response (e.g., TLRs) and the properties of complex decreased in enterocytes (ENR largely) and increased in mixtures of secreted components, including multiple PCs and EECs, both of which occur predominantly in the antimicrobial proteins. Wnt-driven condition ENR + CD (CI-PCs). In total, this suggests that Wnt is necessary for ISCs to commit to PC Conclusions and EEC lineages and that future experimentation with The generation of comprehensive cellular atlases from specific synthetic Wnt ligands [57] may prove fruitful in humans and model organisms will revolutionize our distinguishing Wnt target genes that discriminatorily yield understanding of complex tissues [3]. Intestinal orga- PCs or EECs. Additionally, it is clear that strong Notch noids have already proven their value in studying human Mead et al. BMC Biology (2018) 16:62 Page 16 of 24 and murine epithelial biology. However, to rigorously Sigma-Aldrich) supplemented with growth factors EGF (E; test hypotheses of basic biological or disease mechanism, 50 ng/mL, Life Technologies), Noggin (N; 100 ng/mL, it will be essential to have simple and reliable protocols PeproTech), and R-spondin 1 (R; 500 ng/mL, PeproTech) for the generation of specialized subsets of cells which and small molecules CHIR99021 (C; 3 μM, LC Laborator- cannot be readily isolated from tissue. The representa- ies) and valproic acid (V; 1 mM, Sigma-Aldrich) was added tiveness of cell states present in organoids and the to each well. ROCK inhibitor Y-27632 (Y; 10 μM, R&D specialized cell types present in vivo [3] is an outstand- Systems) was added for the first 2 days of culture. Cells ing question with implications in mucosal immunology, were cultured at 37 °C with 5% CO , and cell culture developmental biology, and translational medicine. Our medium was changed every other day. After 6 days of single-cell genomics approach provides compelling evi- culture, crypt organoids were isolated from Matrigel by dence that organoid-derived cell populations must be mechanical dissociation. Isolated organoids were resus- validated to ensure physiological relevance, and add- pended in TrypLE Express (Life Tech) to dissociate into sin- itionally provides a rational framework for identifying gle cells, then replated in Matrigel with ENR + CV + Y cell states and their potential upstream drivers to modu- media for 2 days. Cells were once again passaged, either late cellular composition. This approach could enable ad- into freezing media (Life Tech) for cryopreservation or vances beyond conventional organoid systems to provide replated at approximately 200 organoids per well (24-well an enriched highly specialized cell population that recapit- plate) for ISC-enriched organoid expansion. ISC-enriched ulates important physiological functions of the intestinal organoids were passaged or differentiated every 6 days in epithelium and could represent an improvement in in the ENR + CV condition. To differentiate, cells were pas- vitro PC culture for the purposes of high-throughput saged as previously described, and crypt culture medium screening, the study of host-microbe interactions, bio- containing growth factors ENR only or ENR + CD (DAPT, engineering (e.g., precision gene editing), and the identifi- 10 μM; Sigma-Aldrich) was added to each well. cation of novel genetic candidates in PC function (e.g., Nupr1). With this framework, we illustrate the power and RNA extraction and qRT-PCR importance of rigorously characterizing the specialized Organoids were isolated from Matrigel in 24-well plates cell types derived in organoids to those defined in ‘atlas-le- following culture as previously described, and pellets vel’ surveys of the intestinal epithelium. were lysed in TRI reagent with RNA extracted according to the manufacturer’s protocol (T9424, Sigma). Resulting Methods RNA pellets were dissolved in UltraPure water and Mice for tissue isolation cDNA synthesis was performed using QuantiTect Re- Proximal small intestine was isolated from C57BL/6 verse Transcription Kit (Qiagen). qPCR reactions were mice of both sexes, aged between 3 and 6 months in all performed using TaqMan Universal Master Mix II (no experiments. UNG), pre-designed TaqMan probes (Additional file 9: Table S5), and 50 ng of sample cDNA (LifeTech). Reac- Bacteria strains tions were performed using an Applied Biosystems Cells were stored at –80 ºC and grown as follows. E. coli 7900HT system. qPCR results were analyzed using RQ strain MG1655 was grown overnight in LB. For experi- manager 1.2 software to obtain CT values used for rela- ments, overnight cultures of MG1655 were resuspended tive quantification to the housekeeping gene Hprt. in M9 supplemented with 0.4% glucose and 0.2% CASa- mino acids. L. lactis strain MG1363 was grown in M17 Confocal imaging of whole cell clusters media supplemented with 0.5% glucose, and E. faecalis ISC-enriched cell clusters (ENR + CV) suspended in 40 μL strain V583 was grown in Brain Heart Infusion media. of Matrigel were seeded onto round coverslips inside a 24-well plate. Cells were treated with ENR + CD, ENR + Crypt culture, enrichment, and differentiation CV, or ENR as previously described. At day 6, organoids Small intestinal crypts were cultured as previously de- were rinsed (PBS0 3X) and fixed on the coverslips by incu- scribed [59]. Briefly, crypts were resuspended in basal bating with 4% paraformaldehyde (PFA) for 30 min at culture medium (Advanced DMEM/F12 with 2 mM room temperature (RT). Gels were blocked and perme- GlutaMAX and 10 mM HEPES; Thermo Fisher Scien- abilized by incubating at RT for 1 hour with 0.1% Triton tific) at a 1:1 ratio with Corning™ Matrigel™ Membrane X-100 and 5% Powerblock in PBS0. Organoids were Matrix – GFR (Fisher Scientific) and plated at the center stained for DEFA and LYZ by incubating with rat of each well of 24-well plates. Following Matrigel anti-mouse Crp1 (Ayabe Lab clone 77-R63, 5 μg/mL, 50X) polymerization, 500 μL of small intestinal crypt culture and rabbit anti-human Lyz (Dako, RRID: AB_2341230, medium (basal media plus 100X N2 supplement, 50X B27 200X) primary antibodies diluted to 10 μg/mL in staining supplement; Life Technologies, 500X N-acetyl-L-cysteine; solution (0.1% Triton X-100 and 10X Powerblock in PBS0) Mead et al. BMC Biology (2018) 16:62 Page 17 of 24 overnight at 4 °C, followed by secondary antibodies Alexa for 20 min to dissociate into single cells. Dissociated Fluor 647 anti-Rabbit IgG (RRID: AB_2563202, 400X) and cells were centrifuged at 300 g for 3 min at 4 °C. The Alexa Fluor 488 anti-Rat IgG (RRID: AB_2563120, 400X) pellet was resuspended in FACS buffer (1% FBS in PBS, diluted in staining solution for 1 h at RT. Actin was stained Thermo Fisher Scientific) and strained into a 5 mL filter with Alexa Fluor 555 Phalloidin (40X) for 20 min, followed cap tube using a 40 μm filter. The cell suspension was by staining of the nucleus with 3 μMDAPIfor 5min. transferred to a flow prep microcentrifuge tube and Coverslips were mounted onto slides with Vectashield and centrifuged at 300 rcf for 3 min. Cell pellets were resus- imaged within 5 days using an Olympus FV2000 confocal pended in a Zombie violet dye (BioLegend, 100X) in microscope. Whole organoid confocal microscopy images FACS buffer for viability staining followed by 1% PFA were processed and analyzed using ImageJ. To determine fixation for 20 min at RT. Pellets were permeabilized thePC puritypercentage, theImageJPoint Picker plugin for 20 min at RT with staining buffer (0.5% Tween-20 wasused to count the numberofnucleito determine the in FACS buffer, Sigma), and co-stained with rabbit total number of cells and to count the number of DEFA- anti-human FITC-Lyz (Dako, RRID: AB_578661, 100X) and LYZ-containing PCs across all z-slices. To investigate and rat anti-mouse APC-CD24 (Biolegend, RRID: cell polarity in whole organoids, individual cells were AB_2565650, 100X) antibodies diluted in staining selected using ImageJ and mean area intensity within buffer for 45 min at RT. Flow cytometry was per- selected cell areas was computed in each z-slice through- formed using a BD LSR II HTS (BD; Koch Institute out the depth of the image across every channel imaged. Flow Cytometry Core at MIT). Initial settings and laser voltages were determined with unstained, single chan- High-resolution single-cell imaging nel stains or secondary-only controls (data not shown). Cell clusters were harvested and rinsed (basal culture Flow cytometry data was analyzed using FlowJo v10.7 media 3X) to remove Matrigel as previously described. software. Briefly, gating was performed as seen in Isolated clusters were resuspended in TrypLE Express Additional file 3: Figure S1F by removing doubles and and incubated at 37 °C for 20 min to dissociate into debris, then selecting the BV421 (viable) cell popula- single cells, then rinsed (basal culture media 2X) and tion; within this population, gating was based on LYZ resuspended in PBS containing magnesium and calcium. and CD24 populations. Pre-coated poly-L-lysine coverslips (Fisher Scientific) were placed into wells of a 24-well plate, a cell suspen- Lysozyme functional secretion assay sion containing approximately 50,000 cells per well was Lysozyme secretion was measured using a Lysozyme added to each well, and the plate was centrifuged at 700 Assay Kit (EnzChek; Thermo Fisher). Briefly, cells sus- rcf for 5 min. PBS supernatant was removed from the pended in Matrigel in 24-well plates were washed (basal wells, and the cells attached to the coverslips were fixed culture media 3X) and either supplemented with 500 μL by incubating with 4% PFA for 30 min at RT. After each of basal culture media or basal culture media plus step, cells were rinsed (PBS 2–5 min 3X). Cells were 10 μM CCh (Sigma Aldrich) for 24 h at 37 °C. Following blocked and permeabilized by incubating at RT for stimulation, culture plates were spun at high speed (> 30 min with permeabilization solution and stained with 2000 g) for 5 min at RT to pellet cell debris and loose DEFA and LYZ by incubating with rat anti-mouse Crp1 Matrigel, and 25 μL of conditioned supernatant was (Ayabe Lab clone 77-R63, 5 μg/mL, 50X) and rabbit removed from the top of each well and quantified as per anti-human Lyz (Dako, RRID: AB_2341230, 200X) pri- the manufacturer’s protocol. mary antibodies diluted in staining solution overnight at 4 °C. Secondary antibodies Alexa Fluor 647 anti-Rabbit Quantification of cell viability, apoptosis, and cytotoxicity IgG (RRID: AB_2563202, 400X) and Alexa Fluor 488 To track proliferation and cell viability, DNA content anti-Rat IgG (RRID: AB_2563120, 400X) diluted in stain- was quantified over the course of differentiation and ing solution were incubated with the coverslips for 1 h CCh stimulation using a CyQUANT Cell Proliferation at RT. Actin was stained with Alexa Fluor 555 Phalloidin Assay Kit (Thermo Fisher) as per the manufacturer’s incubated for 20 min at RT, and the nucleus was stained protocol. Briefly, culture media was aspirated from each with DAPI by incubating at RT for 5 mins. Coverslips well, and the wells were washed (PBS 3X). Gels were were mounted on to slides with Vectashield and imaged then mechanically dissociated into PBS, the contents within 48 h using an Applied Precision DeltaVision transferred into a Falcon tube, centrifuged at 300 rcf for Microscope. 3 min at 4 °C, and the pellet resuspended in PBS to wash. Tubes were centrifuged at 300 rcf for 5 min at Flow cytometry 4 °C, and the pellet was resuspended in 1 mL of assay Cell clusters were isolated from Matrigel as previously working solution (20X cell-lysis buffer, 400X GR dye in described and resuspended in TrypLE Express at 37 °C DI water); 200 μL of samples and DNA standards were Mead et al. BMC Biology (2018) 16:62 Page 18 of 24 plated in triplicate in a black 96-well plate, shaken for concentration of 1.0% and subsequently desalted on 5 min, then fluorescence was measured on a plate 10 mg OASIS HLB solid phase columns (Waters). reader (480 nm/520 nm). From each condition (n = 8), 50 μg aliquots of the Ng For ENR/ENR + CD co-culture, ISC-enriched orga- KD dried tryptic peptides were reconstituted in 100 mM noids (ENR + CV) were differentiated in ENR and ENR HEPES at pH 8.0 to a final concentration of 1.0 mg/mL. + CD and isolated as previously described. The cell The peptides were labeled with TMT-10 isobaric mass pellets were counted and resuspended in basal culture tag reagent according to the manufacturer’s instructions medium, mixed at 0:100, 25:75, 50:50, 75:25, and 100:0% (ThermoFisher Scientific). The peptides were labeled at ENR:ENR + CD ratios (number of clusters), and plated a 1:8 ratio of peptide to TMT reagent, followed by 1 h as previously described in Matrigel in a 96-well plate at incubation at RT with bench-top shaking at 850 rpm. approximately 50 clusters per well in ENR media. After After incubation, a 1.0 μg aliquot of labeled tryptic 24 and 48 h of co-culture, viability versus cytotoxicity peptide was removed from each labeled condition, and caspase activation were assessed using ApoTox-Glo desalted with C18 stage tips, and analyzed via LC-MS/ Triplex Assay (Promega) according to the manufac- MS using a Thermo Fisher Q Exactive Plus Hybrid Mass turer’s protocol. Briefly, 20 μLof ‘V/C reagent’ (10 μL Spectrometer (QE-Plus) coupled to a Thermo Fisher each of GF-AFC and bis-AAF-R110 substrates in 2.0 mL EASY-nLC 1000 liquid chromatograph to ensure iso- of assay buffer) were added to all wells and mixed by baric label incorporation ≥ 95%. An additional 1.0 μgof orbital shaking at 500 rpm for 30 s. After 30 min of labeled tryptic peptide was removed from each channel, incubation at 37 °C, fluorescence was measured on a mixed together, desalted on a C18 stage tip, and ana- plate reader (400 nm/505 nm for viability and 485 nm/ lyzed via LC-MS to ensure equal relative protein loads. 520 nm for cytotoxicity). Caspase-Glo 3/7 reagent (100 μL) During these quality control steps, the labeled peptides was then added to all wells and mixed by orbital shaking at were stored, unquenched, at –80 °C. After validation, 500 rpm for 30 s. After 30 min of incubation at RT, lumi- each channel was quenched with a 5% hydroxylamine nescence was measured on a plate reader. solution to a final sample concentration of 0.3% to quench any unbound isobaric tags. The corresponding eight channels were mixed together for a total amount Bacteria co-culture of 400 μg of labeled tryptic peptides. The labeled peptide For bacteria co-culture, ISC-enriched cells (ENR + CV) mixture was dried down in a speedvac and subsequently were differentiated in ENR and ENR + CD as previously desalted on 30 mg of OASIS HLB solid phase column described. After 6 days of differentiation, cell clusters (Waters). were isolated as previously described. The cell pellet was The dried, labeled peptides were fractionated into 24 resuspended in basal culture medium and plated in sus- fractions by basic reversed-phase using an Agilent pension in a 96-well plate at approximately 150 clusters Zorbax 300 Å, 4.6 mm × 250 mm Extend-C18 column per well. A 1:1 volume of bacteria in respective media on an Agilent 1100 Series HPLC instrument (Agilent (see ‘Bacterial strains’, above; in exponential growth, as Technologies) to decrease sample complexity and in- confirmed by plate reader OD) was added, and bacterial crease the dynamic range of detection. Solvent A (2% growth was measured by serial plating (CFU) after a 4 h acetonitrile, 5 mM ammonium formate, pH 10), and a incubation. Results for bacteria co-culture were normal- non-linear increasing concentration of solvent B (90% ized to no cell (bacteria only) controls. acetonitrile, 5 mM ammonium formate, pH 10) was used as the mobile phase with a flow rate of 1 mL/min Mass spectrometry proteomics sample preparation, through the column. A non-linear gradient with increas- sequencing, and quantification ing percentages of solvent B with four different slopes Organoid cell pellets were isolated from Matrigel with was used (0% for 7 min; 0% to 16% in 6 min; 16% to mechanical dissociation and washed (cold PBS 5X) to 40% in 60 min; 40% to 44% in 4 min; 44% to 60% in remove residual extracellular protein. Proteins were ex- 5 min; 60% for 14 min), and the eluted peptides were tracted from cell pellets with 8 M urea (Sigma), reduced collected in a Whatman polypropylene 2 mL 96-well with 5 mM DTT (Thermo Fisher Pierce) for 45 min, plate (Whatman). The 96 fractions were concatenated alkylated with 10 mM IAA (Sigma) for 45 min in the down to 25 fractions. dark, and double digested with both Lysyl Endopeptidase The global proteome (25 fractions) was analyzed by ‘LysC’ (Wako) and trypsin (Promega) overnight at RT. A LC-MS/MS using the same system described above. –1 small aliquot of cellular lysate was removed from each Peptides were separated at a flow rate of 200 nL min sample for protein quantification via the Pierce™ BCA on a capillary column (Picofrit with a 10-μm tip opening Protein Assay Kit (Pierce). After proteolytic digestion, and 75 μm diameter, New Objective, PF360-75-10-N-5) the samples were quenched using formic acid to a final packed at the Broad Institute with 20 cm of C18 1.9 μm Mead et al. BMC Biology (2018) 16:62 Page 19 of 24 silica beads (1.9 μm ReproSil-Pur C18-AQ medium, Dr. Therefore, in this step of normalization, we normalized Maisch GmbH, r119.aq). Injected peptides were sepa- the median log2 ratio for each ratio column so that the –1 rated at a flow rate of 200 nL min with a linear 84 min median log2 ratio was zero. To robustly and confidently gradient from 100% solvent A (3% acetonitrile, 0.1% for- detect real differential peptides and proteins in the mic acid) to 30% solvent B (90% acetonitrile, 0.1% formic TMT-labeled experiment, we performed a moderated t acid), followed by a linear 9 min gradient from 30% solv- test [63, 64]. Unlike the standard t test, which is not ent A to 90% solvent B for a total of 110 min. The robust for small numbers of samples, the moderated t test QE-Plus instrument was operated in the data-dependent uses an empirical Bayes approach that ‘moderates’ vari- mode acquiring higher-energy collisional dissociation ance estimates for peptides (i.e., shrunk towards a com- tandem mass spectrometry (HCD MS/MS) scans (Reso- mon value), thereby significantly improving the stability of lution = 35,000) for TMT-10 on the 12 most abundant variance estimates for individual peptides. The p values re- ions using an MS1 ion target of 3 × 10 ions and an MS2 ported by the moderated t test were adjusted for multiple target of 5 × 10 ions. The maximum ion time used for testing using the Benjamini–Hochberg FDR method [64]. the MS/MS scans was 120 ms; the HCD-normalized collision energy was set to 31; the dynamic exclusion Proteome pathway and network analysis time was set to 20 s, and the peptide-match preferred Using the identified and quantified proteins from the setting was enabled. TMT-10 labeling experiment, multiple pathway and net- work analyses were performed. Sample correlations were Protein and peptide identification and quantification represented as r-values and determined using GraphPad Peptide spectrum matching and protein identification was Prism version 7.0a. To elucidate potential transcriptional performed using Agilent Technologies SM software pack- drivers of proteome structure, we performed GSEA age (developed at the Broad Institute). In SM, FDRs are v3.0b2 [37, 38] using the full rank-ordered proteome calculated at three different levels, namely spectrum, dis- against the transcription factor target gene set database tinct peptide, and distinct protein. Peptide FDRs are calcu- (v5.2 MSigDB) [39], then performed enrichment map lated in SM using essentially the same pseudo-reversal visualization using GSEA-P-based implementation and strategy evaluated by Elias and Gygi [60] and shown to Cytoscape v3.4.0 [40, 41], with a moderately conservative perform the same as library concatenation. A false distinct cutoff (p < 0.005 and FDR < 0.075) and an overlap coeffi- protein ID occurs when all the distinct peptides that group cient of 0.2. To assess the functional and compartmen- together to constitute a distinct protein have a deltaFor- tal functions associated with the ENR + CD-enriched wardReverseScore ≤ 0. We adjusted the settings to provide proteome and ENR-enriched proteome, we used DA- peptide FDR of 1–2% and protein FDR of 0–1%. SM also VID v6.8 [65, 66] and the gene ontology database, look- carries out sophisticated protein grouping using the ing only at experimentally verified associations within methods previously described [61]. Only proteins with biological processes, cellular compartments, and mo- more than two peptides and at least two TMT ratios in lecular function against a background set of all 8015 each replicate were counted as being identified and quan- quantified proteins. tified. Additionally, we added the capability to flag poten- tially unreliable TMT quantification results based on Single-cell RNA-sequencing detection of more than one precursor in the selection win- A single-cell suspension was obtained from organoids cul- dow for MS/MS. The precursor ion flagging was similar tured under ENR + CV, ENR, and ENR + CD conditions to that recently reported [62], but was carried out post for 6 days as described above. We utilized the Seq-Well data acquisition. As an output, SM generates protein and platform for massively parallel scRNA-seq to capture tran- peptide reports for downstream differential regulation, scriptomes of single cells on barcoded mRNA capture pathway, and network analysis. Prior to comprehensive beads. Full methods on implementation of this platform differential marker, pathway, and network analysis with are available in Gierahn et al. [28]. In brief, 20,000 cells the SM-generated protein reports, we ensured that the from one organoid condition were loaded onto one array data was of high quality and has been properly normal- containing 100,000 barcoded mRNA capture beads. The ized. The first level of normalization was accomplished by loaded arrays containing cells and beads were then sealed guaranteeing that an equivalent amount of peptide (50 μg using a polycarbonate membrane with a pore size of per) was labeled for each of the 10 TMT channels. Once 0.01 μm, which allows for exchange of buffers but retains the SM reports were generated, we calculated the median biological molecules confined within each microwell. ratios for each of the channels where the denominator of Subsequent exchange of buffers allows for cell lysis, tran- the ratio was a predetermined TMT channel signifying the script hybridization, and bead recovery before performing control condition. The underlying assumption was that reverse transcription en masse. Following reverse tran- the null distribution is centered at zero in log2 space. scription and exonuclease treatment to remove excess Mead et al. BMC Biology (2018) 16:62 Page 20 of 24 primers, PCR amplification was carried out using KAPA total number of ENR + CV, ENR, and ENR + CD cells in- HiFi PCR Mastermix with 2000 beads per 50 μLreaction cluded in the analysis was 985, 2544, and 2382, respect- volume. Six libraries (totaling 12,000 beads) were then ively, with quality metrics for nGene, nUMI, and pooled and purified using Agencourt AMPure XP beads percentage of ribosomal and mitochondrial genes re- (Beckman Coulter, A63881) by a 0.6X SPRI followed by a ported in Additional file 8: Figure S4. We then per- 0.7X SPRI and quantified using Qubit hsDNA Assay formed principal component analysis over the list of (Thermo Fisher). Libraries were constructed using the variable genes. For both clustering and t-stochastic Nextera Tagmentation method on a total of 800 pg of neighbor embedding (tSNE), we utilized the first 12 pooled cDNA library from 12,000 recovered beads. Tag- principal components based on the elbow method, as mented and amplified sequences were purified at a 0.6X upon visual inspection of genes contained within, each SPRI ratio yielding library sizes with an average distribu- contributing to important biological processes of intes- tion of 650–750 base pairs in length as determined using tinal cells. We used FindClusters with a resolution of the Agilent hsD1000 Screen Tape System (Agilent Gen- 1.35 and 1000 iterations of tSNE to identify 14 clusters omics). Arrays were sequenced with an Illumina 75 Cycle across the three input samples. To identify genes NextSeq500/550v2 kit at a final concentration of 2.8 pM. which defined each cluster, we performed a ROC test The read structure was paired end with Read 1 starting implemented in Seurat with a threshold set to an from a custom read 1 primer containing 20 bases with AUC of 0.60. a 12 bp cell barcode and 8 bp UMI and Read 2 being 50 bases containing transcript information (Additional Transcriptional scoring file 10: Table S6 for primers used). To determine the fractional contribution to a cell’s tran- scriptome of a gene list, we summed the total log(scaled Single-cell RNA-sequencing computational pipelines and UMI + 1) expression values for genes within a list of analysis interest and divided by the total amount of scaled UMI Read alignment was performed as in Macosko et al. [67]. detected in that cell giving a proportion of a cell’s tran- Briefly, for each NextSeq sequencing run, raw sequen- scriptome dedicated to producing those genes. From our cing data was converted to demultiplexed FASTQ files proteomic screen, we took a list of upregulated proteins using bcl2fastq2 based on Nextera N700 indices corre- (249) or downregulated proteins (212) that were de- sponding to individual samples/arrays. Reads were then tected within our single-cell RNA-sequencing data. To aligned to mm10 genome using the Galaxy portal main- determine the relationship to in vivo PCs and EECs, we tained by the Broad Institute for Drop-Seq alignment took reference data from two Seq-Well experiments run using standard settings. Individual reads were tagged on epithelial cells dissociated from the ileal region of the according to the 12 bp barcode sequencing and the 8 bp small intestine of two C57BL/6 J mice run in separate UMI contained in Read 1 of each fragment. Following experiments. The ileum was first rinsed in 30 mL of ice alignment, reads were binned onto 12 bp cell barcodes cold PBS and allowed to settle. The segment was then and collapsed by their 8 bp UMI. Digital gene expression sliced with scissors and transferred to 10 mL of epithe- matrices (e.g., cell by gene tables) for each sample were lial cell solution (HBSS Ca/Mg-Free 10 mM EDTA, obtained from quality filtered and mapped reads and 100 U/mL penicillin, 100 μg/mL streptomycin, 10 mM UMI-collapsed data, were deposited in GSE100274, and HEPES, 2% FCS (ThermoFisher)) freshly supplemented were utilized as input into Seurat (https://github.com/ with 200 μL of 0.5 M EDTA. The epithelial separation satijalab/seurat) for further analysis [68]. from the underlying lamina propria was performed for To analyze ENR + CV, ENR, and ENR + CD organoids 15 min at 37 °C in a rotisserie rack with end-over-end together, we merged UMI matrices across all genes de- rotation. The tube was then removed and placed on ice tected in any condition and generated a matrix retaining immediately for 10 min before shaking vigorously 15 all cells with at least 1000 UMI detected. This table was times. Visual macroscopic inspection of the tube at this then utilized to setup the Seurat object in which any cell point yielded visible epithelial sheets, and microscopic with at least 400 unique genes was retained and any examination confirmed the presence of single-layer gene expressed in at least five cells was retained. The sheets and crypt-villus structures. The epithelial fraction object was initiated with log-normalization, scaling, and was spun down at 400 g for 7 min and resuspended in centering set to True. Before performing dimensionality 1 mL of epithelial cell solution before transferring to a reduction, data was subset to include cells with less than 1.5 mL Eppendorf tube to minimize the time spent 8000 UMI, and a list of 1676 most variable genes was centrifuging. Cells were spun down at 800 g for 2 min generated by including genes with an average normal- and resuspended in TrypLE Express for 5 min in a 37 °C ized and scaled expression value greater than 0.14 and bath followed by gentle trituration with a P1000 pipette. with a dispersion (variance/mean) greater than 0.4. The Cells were spun down at 800 g for 2 min and Mead et al. BMC Biology (2018) 16:62 Page 21 of 24 resuspended in ACK lysis buffer (ThermoFisher) for Additional file 7: Figure S3 and Additional file 8:FigureS4; 3 min on ice to remove red blood cells and dying cells. n = 8 single-well replicates from one and five biological Cells were spun down at 800 g for 2 min and resus- donors for data reported in Fig. 5a, b and c, respectively; pended in 1 mL of epithelial cell solution and placed on n = 13 co-culture well replicates randomly selected with- ice for 3 min before triturating with a P1000 pipette and out replacement from four donors for data reported in filtering into a new Eppendorf through a 40 μmcell Fig. 5d; n = 6 well replicates (two per three biological strainer (Falcon/VWR). Cells were spun down at 800 g for donors) in Fig. 5e; n = 3 biological donors in Fig. 6c,and 2 min and then resupended in 200 μL of epithelial cell so- n = 5 biological donors in Fig. 6d. For scores in single-cell lution and placed on ice for counting. Single-cell RNA-seq data, we report effect sizes in addition to statistical signifi- data was then generated as described in the ‘Single-cell cance as an additional metric for the magnitude of the RNA-sequencing’ and ‘Single-cell RNA-sequencing com- effect observed. The calculation was performed as Cohen’s putational pipelines and analysis’ sections in Methods.To d, where effect size d =(Mean -Mean )/(SD pooled). 1 2 generate PC and EEC signatures, we ran unbiased SNN-graph based clustering, performed a ROC test, identi- Additional files fied the two mature PC and EEC clusters, and reported all genes with an AUC above 0.60, using all genes with an Additional file 1: Table S1. Derived gene list of the top defining genes from in vivo ileal small intestine PCs and EECs captured on the Seq-Well AUC above 0.65 for scoring within each cluster (gene lists platform. (XLSX 355 kb) in Additional file 1: Table S1), representing any gene with Additional file 2: Table S2. Reference gene lists used in single-cell enrichment in PCs and EECs. These lists capture genes analyses. (XLSX 24 kb) which are enriched in PC (Lyz-high) and EECs (Chga-high) Additional file 3: Figure S1. Image analysis of cell clusters and flow and separate them from the rest of the cells present in in- cytometry. A Bright-field microscopy after 6 days of ENR + CD culture shows annular morphology and darkened lumen of cell clusters consistent with testinal epithelium. For pathway analysis, we inspected cu- presence of granule-rich cells. B Percentage of total cells that are LYZ+ and rated gene lists deposited in the GSEA platform and used DEFA+ following 6 days of ENR, ENR + CV, and ENR + CD culture (from cell KEGG-derived Wnt and Reactome-derived Notch and Re- counting of whole clusters) (n = 3 minimum biological replicates, one-way ANOVA with multiple comparison test versus ENR, **** adj. p < 0.0001). spiratory Electron Transport Chain signatures (Additional C Collapsed z-stack of whole cluster with individual cells highlighted (1–3) file 2: Table S2). In vivo transcription factors for PCs and following 6 days of ENR + CD, stained for LYZ and DEFA and counterstained EECs were determined by matching the PC and EEC signa- with DAPI and for actin (phalloidin). (1–3) Normalized mean-area intensity versus z-axis depth profiles of representative individual LYZ+/DEFA+ ture gene sets with transcription factors from the Riken co-staining cells. D Representative flow cytometry of ENR and ENR + CD Transcription Factor Database (TFdb - http://genome.g- at 6 days with distinct populations of CD24+ and LYZ+ cells indicative sc.riken.jp/TFdb/), and then including only those TFs of phenotypic PCs. E Representative gating for flow cytometry, including removal of doublets and non-viable cells in final gating. F Percentage of which were robustly identified in the proteome dataset. viable cells (membrane impermeable) over time of ENR versus ENR + CD culture. (PDF 3217 kb) Quantification and statistical analysis Additional file 4: Figure S2. Proteomic pipeline, sample-to-sample All statistical analyses were performed using GraphPad comparison, and insights from the in vitro PC proteome. A Schematic of proteomic analysis for samples: culture, collection, lysis, reduction and Prism v7.0a, Seurat implemented in RStudio, and Agilent alkylation, proteolytic digestion, labeling of peptides with isobaric mass tag Technologies Spectrum Mill software package. All graphs reagents (Tandem Mass Tags, TMT10-plex; Thermo), off-line fractionation by with n > 6 show mean ± SEM, unless otherwise noted, basic reverse phase chromatography, analysis of fractions by LC-MS/MS, identification of peptides and proteins using Spectrum Mill software whereas graphs with n < 6 show mean and individual repli- (Agilent), and statistical analysis of the resulting data (moderated t test) to cate values. Unpaired two-tail t test and two-way ANOVA identify confidently differential proteins. B Proteome sample correlation with Dunnett’s multiple comparison test (reported as adj. p between all biological (n = 2) and technical (n = 2/biological) replicates. C ENR + CD-enriched proteins are well-annotated in the gene ontology value) were used to assess statistical significance as appro- (GO) database and show robust enrichment for functions and compartments priate and, unless otherwise noted, * indicates p <0.05, ** p of secretory cells determined by fold enrichment vs. FDR using DAVID. D ENR- < 0.01 *** p < 0.001, **** p < 0.0001, and ns non-significant. enriched proteins are well annotated in the GO database and show enrichment for functions and compartments of transcriptionally and In each experiment, tissues were isolated from multiple translationally active cells determined by fold enrichment vs. FDR mice housed in the same facility with each mouse providing using DAVID. e GSEA enrichment map of transcription factors linked tissue designated as a distinct biological donor: n =3 to ENR + CD- and ENR-enriched proteins following a moderately conservative cutoff of p < 0.005, FDR < 0.075, and overlap coefficient donor-averaged values of four technical replicates for data of 0.2. (PDF 483 kb) reported in Fig. 2b; n = 3 donor of two technical replicates Additional file 5: Table S3. Detected and quantified in vitro Proteome. for data reported in Fig. 2e (with exception of day 8, n =2 (XLSX 919 kb) donor-averaged values of two technical replicates) and Additional file 6: Table S4. Complete list of DAVID enrichments. Additional file 3:FigureS1F; n = 4 (two technical replicates (XLSX 51 kb) from two biological donors each) for data reported in Additional file 7: Figure S3. Quality metrics for single-cell RNA sequencing. A Total gene number of cells maintained in analyses with a lower cutoff of Figs. 2f–h, and Additional file 4:FigureS2; n = 1 biological n = 400 unique genes per cell. Total unique molecular identifiers (UMIs) used donor for in vitro data reported in Figs. 3 and 4 and Mead et al. BMC Biology (2018) 16:62 Page 22 of 24 Competing interests as the basis for cell-by-gene tables collapsed to UMI as input into Seurat with JMK, RSL, and XY hold equity in Frequency Therapeutics, a company that has the lower bound representing n = 400 unique genes and the upper bound an option to license IP generated by JMK, RSL, and XY and that may benefit representing 8000 UMIs. Note: Clusters ENR + CV-3, ENR + CV-4, and ENR-1 had financially if the IP is licensed and further validated. The interests of JMK, RSL, significantly higher levels of genes and UMIs and, intriguingly, were also the and XY were reviewed and are subject to a management plan overseen by three clusters with the highest levels of Lgr5 (see Fig. 5a), indicating that stem their institutions in accordance with their conflict of interest policies. cells may contain larger contents of RNA, as they are in a biosynthetic state before differentiation and maturation. B Violin plot of expression contribution to a cell’s transcriptome of mitochondrial and ribosomal Publisher’sNote genes across identified subsets. (PDF 1153 kb) Springer Nature remains neutral with regard to jurisdictional claims in Additional file 8: Figure S4. Signaling pathways and processes associated published maps and institutional affiliations. with in vitro PC enrichment (Additional file 9: Table S5 for reference gene lists). A Violin plot of expression contribution to a cell’s transcriptome of Wnt Author details pathway genes (Additional file 2: Table S2; activated by CHIR99021) across 1 Division of Engineering in Medicine, Department of Medicine, Brigham & clusters as percent of transcriptome. B Violin plot of expression contribution 2 Women’s Hospital, Harvard Medical School, Boston, MA, USA. Harvard-MIT to a cell’s transcriptome of Notch pathway genes (Additional file 2:Table S2; 3 Division of Health Sciences & Technology, Cambridge, MA, USA. Koch inhibited by DAPT) across clusters as percent of transcriptome. C Violin plot Institute for Integrative Cancer Research, MIT, Cambridge, MA, USA. Harvard of expression contribution to a cell’s transcriptome of respiratory electron 5 Stem Cell Institute, Cambridge, MA, USA. Broad Institute of Harvard and MIT, transport gene set (Additional file 2:Table S2) acrossclustersaspercent of 6 Cambridge, MA, USA. Institute for Medical Engineering and Science, MIT, transcriptome. (PDF 888 kb) 7 Cambridge, MA, USA. Department of Chemistry, MIT, Cambridge, MA, USA. 8 9 Additional file 9 : Table S5. TaqMan gene expression assays used for Department of Chemical Engineering, MIT, Cambridge, MA,, USA. Ragon qRT-PCR. (DOCX 13 kb) Institute of MGH, MIT and Harvard, Cambridge, MA, USA. Divisions of Infectious Diseases and Gastroenterology, Boston Children’s Hospital and Additional file 10: Table S6. SeqWell reverse transcription and library Harvard Medical School, Boston, MA, USA. Wyss Institute for Biologically preparation primers. (DOCX 13 kb) Inspired Engineering, Harvard University, Boston, MA, USA. Department of Biological Engineering, MIT, Cambridge, MA, USA. Synthetic Biology Center, MIT, Cambridge, MA, USA. Center for Microbiome Informatics and Acknowledgements Therapeutics, MIT, Cambridge, MA, USA. We thank the Koch Institute of MIT Core facilities for assisting with microscopy, flow cytometry, and pilot mass spectrometry experiments, and Received: 30 March 2018 Accepted: 8 May 2018 for the partial support from the Cancer Center Support (core) Grant P30- CA14051 from the NCI. We also thank R Blumberg, N Krupka, J Matute, and Y Shao for critically reviewing the manuscript. Furthermore, we appreciate the contribution of anti-crypitidin (DEFA) antibodies from Tokiyoshi Ayabe References (Hokkaido University). 1. Clevers H. Modeling development and disease with organoids. Cell. 2016; 165:1586–97. https://doi.org/10.1016/j.cell.2016.05.082 Funding 2. Prakadan SM, Shalek AK, Weitz DA. Scaling by shrinking: empowering This research was supported by both an Innovator award and Breakthrough single-cell “omics” with microfluidic devices. Nat Rev Genet. 2017;18:345–61. award to JMK from the Kenneth Rainin Foundation, the US National Institutes of 3. Haber AL, Biton M, Rogel N, Herbst RH, Shekhar K, Smillie C, et al. A single- Health (NIH) grant DE013023 to RL, NIH grant HL095722 to JMK, and a ‘shark tank’ cell survey of the small intestinal epithelium. Nature. 2017;551(7680):333–9. pilot grant from Brigham and Women’sHospitaltoJMK.BEM wassupported by 4. Grün D, Lyubimova A, Kester L, Wiebrands K, Basak O, Sasaki N, et al. Single- the Graduate Research Fellowship Program from the US National Science cell messenger RNA sequencing reveals rare intestinal cell types. Nature. Foundation (NSF). AKS was supported by the Searle Scholars Program, the 2015;525:251–5. Beckman Young. 5. The HCA Consortium. The human cell atlas white paper. 2017; https://www. Investigator Program, a Sloan Research Fellowship in Chemistry, NIHgrants humancellatlas.org/files/HCA_WhitePaper_18Oct2017.pdf. Accessed 29 May 2018. 1DP2OD020839, 2U19AI089992, 2R01HL095791, 1U54CA217377,2P01AI039671, 6. Tanay A, Regev A. Scaling single-cell genomics from phenomenology to 5U24AI118672, 2RM1HG006193, 1R33CA202820, 1R01HL126554,1R01DA046277 mechanism. Nature. 2017;541:331–8. and 1R01AI138546, and Bill and Melinda Gates Foundationgrants OPP1139972, 7. Satija R, Shalek AK. Heterogeneity in immune responses: From populations OPP1137006, and OPP1116944. JOM is a Damon Runyon Fellow supported by to single cells. Trends Immunol. 2014;35:219–29. https://doi.org/10.1016/j.it. the HHMI and Damon Runyon Cancer Research Foundation DRG-2274-16. 2014.03.004 8. Lancaster MA, Knoblich JA. Organogenesis in a dish: Modeling development and disease using organoid technologies. Science. 2014;345:1247125. Availability of data and materials 9. Schwank G, Koo BK, Sasselli V, Dekkers JF, Heo I, Demircan T, et al. The accession number for the proteomics data reported in this paper is Functional repair of CFTR by CRISPR/Cas9 in intestinal stem cell organoids MassIVE: MSV000081369. The accession number for the single-cell RNA of cystic fibrosis patients. Cell Stem Cell. 2013;13:653–8. https://doi.org/10. sequencing data reported in this paper is GEO: GSE100274. 1016/j.stem.2013.11.002 10. Drost J, van Boxtel R, Blokzijl F, Mizutani T, Sasaki N, Sasselli V, et al. Use of Authors’ contributions CRISPR-modified human stem cell organoids to study the origin of Conceptualization, BEM, JOM, XY, AKS and JMK; Methodology, BEM, XY, JMK, mutational signatures in cancer. Science. 2017;238:eaao3130. JOM, AKS, RA; Formal Analysis, BEM, JOM, MJS, MAM, Investigation, BEM, APB, 11. Molodecky NA, Soon IS, Rabi DM, Ghali WA, Ferris M, Chernoff G, et al. LEL, JOM, MJS, DAA, PB, TKH, MHW, SRN; Resources, SRN, RL, JMK, SAC, AKS, Increasing incidence and prevalence of the inflammatory bowel diseases Writing – Original Draft, BEM, LEL, JOM, APB; Writing – Review & Editing, with time, based on systematic review. Gastroenterology. 2012;142:46–54. BEM, JMK, LEL, JOM, PB, XY, RA, SAC, JJC, AKS, RL; Visualization, BEM, JOM, e42; quiz e30 LEL, APB; Supervision, JMK, BEM, JJC, AKS, SAC; Project Administration, JMK, 12. Wehkamp J, Salzman NH, Porter E, Nuding S, Weichenthal M, Petras RE, BEM; Funding Acquisition, JMK, AKS, RL, BEM. All authors read and approved et al. Reduced Paneth cell α-defensins in ileal Crohn’s disease. Proc Natl the final manuscript. Acad Sci U S A. 2005;102:18129–34. 13. Ireland H, Houghton C, Howard L, Winton DJ. Cellular inheritance of a Ethics approval and consent to participate Cre-activated reporter gene to determine Paneth cell longevity in the Isolation of tissues, housing, and maintenance of animal colonies were all murine small intestine. Dev Dyn. 2005;233:1332–6. performed under the approval and guidelines of the MIT Committee on 14. Sato T, van Es JH, Snippert HJ, Stange DE, Vries RG, van den Born M, et al. Animal Care or the Harvard Digestive Disease Center and Brigham and Paneth cells constitute the niche for Lgr5 stem cells in intestinal crypts. Women’s Gnotobiotic Core. Nature. 2011;469:415–8. Mead et al. BMC Biology (2018) 16:62 Page 23 of 24 15. Clevers HC, Bevins CL. Paneth cells: maestros of the small intestinal crypts. 40. Merico D, Isserlin R, Stueker O, Emili A, Bader GD. Enrichment map: a Annu Rev Physiol. 2013;75:289–311. network-based method for gene-set enrichment visualization and 16. Xavier RJ, Podolsky DK. Unravelling the pathogenesis of inflammatory bowel interpretation. PLoS One [Internet]. 2010;5:e13984. Available from: disease. Nature. 2007;448:427–34. http://www.ncbi.nlm.nih.gov/pubmed/21085593 17. Khor B, Gardet A, Xavier RJ. Genetics and pathogenesis of inflammatory 41. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. bowel disease. Nature. 2011;474:307–17. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504. 18. Liu T-C, Gurram B, Baldridge MT, Head R, Lam V, Luo C, et al. Paneth cell 42. Stringari C, Edwards RA, Pate KT, Waterman ML, Donovan PJ, Gratton E. defects in Crohn’s disease patients promote dysbiosis. JCI Insight. 2016;1:1–15. 19. Adolph TE, Tomczak MF, Niederreiter L, Ko H-J, Böck J, Martinez-Naves E, Metabolic trajectory of cellular differentiation in small intestine by Phasor et al. Paneth cells as a site of origin for intestinal inflammation. Nature. Fluorescence Lifetime Microscopy of NADH. Sci Rep. 2012;2:568. 2013;503:272–6. 43. Rodríguez-Colman MJ, Schewe M, Meerlo M, Stigter E, Gerrits J, Pras-Raves 20. Kobayashi KS, Chamaillard M, Ogura Y, Henegariu O, Inohara N, Nuñez G, M, et al. Interplay between metabolic identities in the intestinal crypt et al. Nod2-dependent regulation of innate and adaptive immunity in the supports stem cell function. Nature. 2017;543(7645):424–7. intestinal tract. Science. 2005;307:731–4. 44. Ayabe T, Satchell DP, Wilson CL, Parks WC, Selsted ME, Ouellette AJ. 21. Kaser A, Blumberg RS. ATG16L1 Crohn’s disease risk stresses the Secretion of microbicidal alpha-defensins by intestinal Paneth cells in endoplasmic reticulum of Paneth cells. Gut. 2014;63(7):1038–9. response to bacteria. Nat Immunol. 2000;1:113–8. 45. Kanamori M, Konno H, Osato N, Kawai J, Hayashizaki Y, Suzuki H. A genome- 22. Kaser A, Lee A-H, Franke A, Glickman JN, Zeissig S, Tilg H, et al. XBP1 links wide and nonredundant mouse transcription factor database. Biochem ER stress to intestinal inflammation and confers genetic risk for human Biophys Res Commun. 2004;322:787–93. inflammatory bowel disease. Cell. 2008;134:743–56. 23. Farin HF, Karthaus WR, Kujala P, Rakhshandehroo M, Schwank G, Vries RGJ, 46. Jia SN, Lin C, Chen DF, Li AQ, Dai L, Zhang L, et al. The transcription factor et al. Paneth cell extrusion and release of antimicrobial products is directly p8 regulates Autophagy in response to palmitic acid stress via a controlled by immune cell-derived IFN-γ. J Exp Med. 2014;211:1393–405. mammalian target of rapamycin (mTOR)-independent signaling pathway. J 24. Wilson SS, Tocchi A, Holly MK, Parks WC, Smith JG. A small intestinal Biol Chem. 2016;291:4462–72. organoid model of non-invasive enteric pathogen-epithelial cell 47. Grasso D, Bintz J, Lomberk G, Molejon MI, Loncle C, Garcia MN, et al. Pivotal interactions. Mucosal Immunol. 2014;8:1–10. role of the chromatin protein Nupr1 in Kras-induced senescence and 25. Foulke-Abel J, In J, Yin J, Zachos NC, Kovbasnjuk O, Estes MK, et al. Human transformation. Sci Rep. 2015;5:17549. enteroids as a model of upper small intestinal ion transport physiology and 48. Cano CE, Hamidi T, Sandi MJ, Iovanna JL. Nupr1: The Swiss-knife of cancer. pathophysiology. Gastroenterology. 2016;150:638–49. e8 J Cell Physiol. 2011;226:1439–43. 49. Imielinski M, Baldassano RN, Griffiths A, Russell RK, Annese V, Dubinsky M, 26. Moon C, VanDussen KL, Miyoshi H, Stappenbeck TS. Development of a primary et al. Common variants at five new loci associated with early-onset mouse intestinal epithelial cell monolayer culture system to evaluate factors inflammatory bowel disease. Nat Genet. 2009;41:1335–40. that modulate IgA transcytosis. Mucosal Immunol. 2013;7:818–28. 27. Basak O, Beumer J, Wiebrands K, Seno H, van Oudenaarden A, Clevers H. 50. Santofimia-Castaño P, Rizzuti B, Pey ÁL, Soubeyran P, Vidal M, Urrutia R, Induced quiescence of Lgr5+ stem cells in intestinal organoids enables et al. Intrinsically disordered chromatin protein NUPR1 binds to the differentiation of hormone-producing enteroendocrine cells. Cell Stem Cell. C-terminal region of Polycomb RING1B. Proc Natl Acad Sci. 2017; 2017;20:177–90. e4 https://doi.org/10.1073/pnas.1619932114. 28. Gierahn TM, Wadsworth MH, Hughes TK, Bryson BD, Butler A, Satija R, et al. 51. Neira JL, Bintz J, Arruebo M, Rizzuti B, Bonacci T, Vega S, et al. Identification Seq-Well: portable, low-cost RNA sequencing of single cells at high of a drug targeting an intrinsically disordered protein involved in pancreatic throughput. Nat Methods. 2017;14:1–8. adenocarcinoma. Sci Rep. 2017;7:39732. 29. Yin X, Farin HF, van Es JH, Clevers H, Langer R, Karp JM. Niche-independent 52. Gulati AS, Shanahan MT, Arthur JC, Grossniklaus E, von Furstenberg RJ, high-purity cultures of Lgr5+ intestinal stem cells and their progeny. Nat Kreuk L, et al. Mouse background strain profoundly influences Paneth Methods. 2014;11:106–12. cell function and intestinal microbial composition. PLoS One. 2012;7: e32403. 30. Yin X, Mead BE, Safaee H, Langer R, Karp JM, Levy O. Engineering stem cell organoids. Cell Stem Cell. 2016;18:25–38. 53. Bevins CL, Salzman NH. Paneth cells, antimicrobial peptides and maintenance of intestinal homeostasis. Nat Rev Microbiol. 31. McLean WJ, Yin X, Lu L, Lenz DR, McLean D, Langer R, et al. Clonal expansion of Lgr5-positive cells from mammalian cochlea and high-purity 2011;9:356–68. generation of sensory hair cells. Cell Rep. 2017;18:1917–29. 54. Zhang Q, Pan Y, Yan R, Zeng B, Wang H, Zhang X, et al. Commensal 32. van Es JH, Jay P, Gregorieff A, van Gijn ME, Jonkheer S, Hatzis P, et al. Wnt bacteria direct selective cargo sorting to promote symbiosis. Nat Immunol. signalling induces maturation of Paneth cells in intestinal crypts. Nat Cell 2015;16(9):918–26. Biol. 2005;7:381–6. 55. Cunliffe RN, Rose FR, Keyte J, Abberley L, Chan WC, Mahida YR. Human 33. VanDussen KL, Carulli AJ, Keeley TM, Patel SR, Puthoff BJ, Magness ST, et al. defensin 5 is stored in precursor form in normal Paneth cells and is Notch signaling modulates proliferation and differentiation of intestinal expressed by some villous epithelial cells and by metaplastic Paneth cells in crypt base columnar stem cells. Development. 2012;139:488–97. the colon in inflammatory bowel disease. Gut. 2001;48:176–85. 56. Beumer J, Clevers H. Regulation and plasticity of intestinal stem cells during 34. Tian H, Biehs B, Chiu C, Siebel CW, Wu Y, Costa M, et al. Opposing activities homeostasis and regeneration. Development. 2016;143:3639–49. of notch and wnt signaling regulate intestinal stem cells and gut homeostasis. Cell Rep. 2015;11:33–42. 57. Janda CY, Dang LT, You C, Chang J, de Lau W, Zhong ZA, et al. Surrogate 35. Buczacki SJ, Zecchini HI, Nicholson AM, Russell R, Vermeulen L, Kemp R, Wnt agonists that phenocopy canonical Wnt and β-catenin signalling. et al. Intestinal label-retaining cells are secretory precursors expressing Lgr5. Nature. 2017;545(7653):234–7. Nature. 2013;495:65–9. 58. Gjorevski N, Sachs N, Manfrin A, Giger S, Bragina ME,Ordóñez-Morán P et al. 36. von Furstenberg RJ, Gulati AS, Baxi A, Doherty JM, Stappenbeck TS, Gracz Designer matrices for intestinal stem cell and organoid culture. Nature 2016; AD, et al. Sorting mouse jejunal epithelial cells with CD24 yields a 539:560–564. https://doi.org/10.1038/nature20168. population with characteristics of intestinal stem cells. AJP Gastrointest Liver 59. Sato T, Vries RG, Snippert HJ, van de Wetering M, Barker N, Stange DE, et al. Physiol. 2011;300:G409–17. Single Lgr5 stem cells build crypt-villus structures in vitro without a 37. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. mesenchymal niche. Nature. 2009;459:262–5. Gene set enrichment analysis: A knowledge-based approach for interpreting 60. Elias JE, Gygi SP. In: Hubbard SJ, Jones AR, editors. Target-Decoy Search genome-wide expression profiles. Proc Natl Acad Sci. 2005;102:15545–50. Strategy for Mass Spectrometry-Based Proteomics. Totoewa: Humana Press; 2010. p. 55–71. http://link.springer.com/10.1007/978-1-60761-444-9. 38. Mootha VK, Lindgren CM, Eriksson K-F, Subramanian A, Sihag S, Lehar J, et al. PGC-1α-responsive genes involved in oxidative phosphorylation are 61. Nesvizhskii AI, Aebersold R. Interpretation of Shotgun proteomic data. Mol coordinately downregulated in human diabetes. Nat Genet. 2003;34:267–73. Cell Proteomics. 2005;4:1419–40. 39. Xie X, Lu J, Kulbokas EJ, Golub TR, Mootha V, Lindblad-Toh K, et al. 62. Phanstiel DH, Brumbaugh J, Wenger CD, Tian S, Probasco MD, Bailey DJ, Systematic discovery of regulatory motifs in human promoters and 3’ UTRs et al. Proteomic and phosphoproteomic comparison of human ES and iPS by comparison of several mammals. Nature. 2005;434:338–45. cells. Nat Methods. 2011;8:821–7. Mead et al. BMC Biology (2018) 16:62 Page 24 of 24 63. Smyth GK. Linear Models and Empirical Bayes Methods for Assessing Differential Expression in Microarray Experiments. Stat Appl Genet Mol Biol. 2004;3:1–26. 64. Benajmini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc B. 1995;57:289–300. 65. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4:44–57. 66. Huang DW, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: Paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009;37:1–13. 67. Macosko EZ, Basu A, Satija R, Nemesh J, Shekhar K, Goldman M, et al. Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets. Cell. 2015;161:1202–14. 68. Satija R, Farrell JA, Gennert D, Schier AF, Regev A. Spatial reconstruction of single-cell gene expression data. Nat Biotechnol. 2015;33:495–502.

Journal

BMC BiologySpringer Journals

Published: Jun 5, 2018

References

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off