Dictionary of human intestinal organoid responses to secreted niche factors at single cell resolution

Domain: biology
Published: 2026-01-09
Editorial outcome: major_revision

Abstract
The intestinal epithelium plays a critical role in health and disease, yet the impact of microenvironmental cues in diseased contexts, such as inflammatory bowel disease (IBD), remains poorly defined. To address this gap, we first benchmarked human colonic organoid injury models against IBD tissue and established a disease-relevant model of inflammation using inflammatory cytokines. Using this system, we built a dictionary of epithelial responses to 79 secreted niche factors at single cell resolution via donor-pooled, multiplexed single cell RNA-sequencing. The comprehensive nature of our atlas allowed us to map relationships between perturbations, infer the function of less characterized ligands, and identify cell type-specific perturbed pathways. Finally, we established the relevance of organoid-derived gene programs by mapping them to single cell and spatial atlases of human IBD tissue. Our resource offers a global view of epithelial responses to microenvironmental cues, offering insights into epithelial homeostasis and repair mechanisms in IBD.

Introduction
The colonic epithelium is a highly regenerative tissue composed of proliferative stem cell-containing crypt domains and differentiated cell types that facilitate absorption and mucus secretion1. During homeostasis, spatially-patterned microenvironmental cues direct stem cell proliferation in the crypts, which supports turnover and differentiation into specialized cell types in the inter-crypt domains2. However, the microenvironment becomes altered in the context of injury or disease, particularly inflammatory bowel disease (IBD), a common intestinal disorder characterized by chronic inflammation3. For instance, stromal populations in the stem cell niche are altered during IBD, which may lead to disrupted signaling to the epithelium4–6.

While much is known about the colonic microenvironment during homeostasis and development1,2,7, it is less clear how niche cues contribute to epithelial repair in a diseased context. First, epithelial regeneration is most commonly studied in models which lack important aspects of disease biology, including mouse models which suffer from species-specific differences8,9, immortalized cell lines or 2D cultures which do not capture the complex 3D architecture of the gut10, or healthy homeostatic organoids that do not model regeneration in the face of injury. Organoids derived from IBD patients have been employed to overcome these issues11, but they often lose their inflammatory phenotype over time in culture12. Second, the lack of systematic interrogation of epithelial cell-type-specific responses to niche cues hinders a global understanding of pathway regulation associated with homeostasis and repair. Most studies to date have focused on characterizing a select few perturbations at a time7,13,14, rendering it difficult to map out divergences and similarities between ligands. Conversely, large phenotypic screens without molecular profiling are also powerful but pose challenges in gaining direct mechanistic understanding15. Moreover, because the intestinal epithelium is heterogenous in nature, understanding cell-type-specific responses at single-cell resolution can provide insight on both mechanism and therapeutic targeting. Indeed, recent scRNA-seq chemical perturbation screens have already begun to demonstrate value in unbiased and systematic perturbation mapping16–18.

To address these gaps, we first developed a physiologically relevant epithelial injury model by characterizing the response of healthy human colon organoids subjected to common epithelial stressors (gamma irradiation, cytokine stimulation, dextran sodium sulfate (DSS), and mechanical injury)19–24 and benchmarking them to human IBD tissue. We discovered an IBD-relevant inflammatory model and used it to systematically profile epithelial responses to 79 secreted niche factors known or predicted to interact with the human colonic epithelium during development, homeostasis, and disease25–27 via scRNA-seq. The breadth of our atlas allowed us to infer cell-type-specific effects of less well-characterized ligands via their relationships to known ligand classes and identify epithelial response gene programs that may be implicated in a return to homeostasis. Overall, we present a valuable resource to probe the diversity and context-specificity of human intestinal epithelial responses to endogenous microenvironmental cues, which can be leveraged for therapeutic target discovery and biomarker selection for IBD.

Results
We developed 4 epithelial injury models starting from homeostatic primary colon organoids derived from 3 human donors. Organoids were cultured in a commercial microwell suspension culture system28 to improve ease of imaging, single-cell dissociation, and homogeneity. We confirmed that microwell-grown organoids contain both stem and differentiated cell types and display transcriptional similarity to both traditional organoids cultured in 3D extracellular matrix domes and to fresh human colon tissue (Supplementary Fig. 1, Supplementary Data 1). Organoids were first differentiated to mimic the complexity of cell types found in vivo and then treated with conventional epithelial stressors including DSS, gamma irradiation, and cytokines (TNFɑ, IFNɣ, IL1β combined - termed cytomix)24,29 across various dosages (Fig. 1A, Supplementary Fig. 2A). We first characterized organoid responses via morphology, degree of cell death, and barrier function, as IBD lesions display increased expression of apoptotic markers and decreased expression of tight junction proteins leading to increased permeability30–32. DSS and cytomix treatment resulted in morphological changes from bright and cystic to dark organoid structures, suggesting stress and cell death (Fig. 1B). While control and irradiated organoids were primarily simple squamous or columnar in architecture, DSS-treated organoids showed epithelial thickening and altered polarity, with Phalloidin staining both the apical and basolateral surfaces. In contrast, cytomix treatment resulted in primarily pseudostratified epithelium, which has been reported for IBD patient-derived organoids33. Out of all treatments, only DSS caused significant decreases in proliferation (Fig. 1B, Supplementary Fig. 2B). While irradiation has previously been shown to decrease proliferation in small intestinal organoids14,21, we did not observe decreases in KI67 expression at 6 Gy or 12 Gy, indicating potential differences in susceptibility between the small and large intestine. Finally, activated Caspase-3/7 staining revealed dose-dependent increases in apoptosis after DSS and cytomix, but not irradiation (Fig. 1C, Supplementary Fig. 2C)34. All injuries except for irradiation led to a significant increase in organoid permeability to 4 kDa FITC dextran compared to untreated organoids, with strong dose-dependent effects in both DSS and cytomix treatment (Fig. 1C, Supplementary Fig. 2C). These results suggest that DSS and cytomix elicit a functional stress response in human colon organoids that may bear relevance to human IBD.

To understand the molecular fidelity of these injury models to human IBD, we performed scRNA-seq after one day of treatment (Fig. 1D). Changes in cell type proportions were observed, including decreased stem and transit-amplifying (TA) cells with DSS, and decreased colonocytes and increased TA/stem cells with cytomix (Fig. 1E). Additionally, a unique cell cluster marked by expression of CXCL8, AREG, and ANXA1 was enriched specifically in organoids recovering from dissociation (Fig. 1F). This population may represent a unique cell state brought on by endogenous organoid formation mechanisms as AREG and ANXA1 have been reported as markers of regeneration and de-differentiation and are associated with the revival stem cell (revSC) signature19,35–38. To determine the disease-relevance of organoid molecular responses, we first conducted a meta-analysis across two published single-cell datasets (colon tissue from healthy, Crohn’s disease, and ulcerative colitis patients) to identify conserved gene programs (GPs) dysregulated between diseased vs. healthy epithelium4,25 (Supplementary Fig. 2D, Supplementary Data 2). Two highly conserved gene programs were found: upregulation of GP IBD1, which is associated with inflammation (DUOX2, DUOXA2, CEACAM6, ANXA1)39 and contains known IBD risk genes40,41, and downregulation of GP IBD2, which is associated with cellular differentiation (CD24, CA1, KRT19) and metabolism (COX8A, NDUFB9, CYCS), particularly oxidative phosphorylation4 and electron transport (Fig. 1G). These changes occur across most epithelial cell types (Supplementary Fig. 2D), suggesting a general heightened immune response and altered metabolic function, as previously described42–47. Additionally, IBD1 and IBD2 appear highly predictive of disease state (Supplementary Fig. 2E) and show similar trends when tested in a third unseen dataset48 (Fig. 1G). We then compared the expression of IBD1 and IBD2 across injury models. IBD1 was upregulated in stem and TA cells under DSS, and even more strongly with cytomix in a dose-dependent manner (Fig. 1H, Supplementary Data 3, 11). Conversely, metabolism-related IBD2 was only downregulated in multiple cell types after cytomix treatment, and to a lesser extent with irradiation. Combined, these results suggest that cytomix stimulation bears the strongest phenotypic and molecular resemblance to human IBD tissue.

We further identified GPs associated with each injury in an unbiased manner (Fig. 2A, Supplementary Fig. 3A, B, Supplementary Data 4, 11). Dissociation upregulated wound healing and apoptosis pathways (GP19: DKK1, ANXA1, ANXA5), processes associated with regeneration. Additionally, dissociation upregulated multiple genes associated with the known revSC signature (CLU, ANXA1, TACSTD2). Cytomix upregulated immune response and interferon pathways (GP16: DUOX2, DUOXA2, CXCL10), irradiation altered DNA damage, apoptosis, and translation pathways (GP8: LY6D, BAX, BCL2L1), while DSS upregulated transcription and cytoskeleton reorganization processes (GP6: THRB, HNF4G, LGR4). We additionally identified a general organoid injury signature (GP2: LDHA, ME1, ENO1) which was downregulated across all injuries and associated with metabolic and biosynthetic processes, and may indicate a broad metabolic defect across injured epithelium (Supplementary Fig. 3C, Supplementary Data 4).

To assess the IBD-relevance of injury-elicited GPs, we evaluated their expression in human IBD tissue. Only cytomix-associated GPs (GP3, GP16) exhibited concordant upregulation in two human IBD datasets4,25 (Supplementary Fig. 4A, Supplementary Data 10). We additionally examined the expression and localization of these GPs in healthy (n=4) vs UC (n=16) patient samples via spatial transcriptomics (CosMx, Supplementary Fig. 4B, Supplementary Data 9). Importantly, cytomix-regulated GP3 was enriched in epithelial cells in UC vs healthy, while most other stressor-induced GPs did not exhibit enrichment in disease (Supplementary Fig. 4C). We also examined genes associated with revSC’s and found that this signature was upregulated in inflamed vs. healthy epithelium, but the revSC signature was weakly expressed in epithelial cells and driven by a few genes via spatial transcriptomic data (Supplementary Fig. 4A, C). Together, this suggests that cytomix treatment best recapitulates IBD-like epithelial gene signatures, while conventional stressors (DSS, irradiation) bear less relevance when applied to human colon organoids.

We further hypothesized that organoid responses to perturbations may be primed by stressors and hence be context-specific. To test this, we performed scRNA-seq on DSS and cytomix-treated organoids perturbed with growth factors with known roles in intestinal epithelial cell maintenance: BMP2, BMP6, EGF, HBEGF, IL22, TGFβ1 (Fig. 2B)13,14,49,50. Notably, we found significantly more differentially expressed genes (DEGs, padj < 0.05, log FC > 1) for BMP6 and BMP2 in cytomix-treated organoids compared to DSS, while DSS-treated organoids exhibited many more DEGs in the presence of EGF and HBEGF (Fig. 2C, Supplementary Fig. 5A, B, Supplementary Data 5). BMP2/6 strongly upregulated canonical BMP target genes (ID1, ID3) and genes involved in metabolic processes (AKR1C2, AKR1C3, HSD11B2) and epithelial cell differentiation (DAB2, KRT20, OVOL2) specifically in the context of cytomix, while EGF ligands elicited pathways involving regulation of cell-cell adhesion and anoikis (CEACAM5, CEACAM6, CAV1) specifically in DSS (Fig. 2D, Supplementary Fig. 5A–C). We hypothesized that these differences could be partially attributed to changes in expression of receptors and inhibitors of both pathways after injury. For example, cytomix treatment upregulated BMP receptor BMPR1A while DSS upregulated EGF-family receptors ERBB2 and ERBB3 and downregulated EGFR-inhibitor LRIG1 (Fig. 2E). Interestingly BMP receptors (BMPR1A, ACVR1) and EGFR are upregulated in IBD vs. healthy controls while ERBB2 is downregulated (Supplementary Fig. 5D). Additionally, response genes may already be activated by the injury itself as cytomix-treated organoids upregulated IL22 response genes while DSS upregulated BMP2 response genes (Supplementary Fig. 5E), thus decreasing the dynamic range to observe perturbation responses. This suggests that the choice of injury model plays an important role in evaluating physiologically relevant epithelial responses, and further strengthens the selection of cytomix as an IBD-relevant model51.

Using the low-dose cytomix model, we generated a molecular dictionary of epithelial responses to secreted niche factors in human colon organoids across 3 donors via multiplexed scRNA-seq (Fig. 3A). After filtering and quality control, we obtained 169,996 cells with 3 donors × 83 perturbations, for a total of 249 independent conditions (Supplementary Fig. 6B). We selected a perturbation space enriched in factors predicted to promote epithelial homeostasis, wound healing, and/or regeneration by performing a meta-analysis of putative ligand-receptor interactions between the gut epithelium and its niche via existing single cell atlases (Supplementary Fig. 6A, Supplementary Data 6)25–27. We selected 70 niche factors from this analysis and 9 additional perturbations with known effects on the intestinal epithelium such as IL2250, DKK152, and ISLR53, for a total of 79 perturbations (36.8% predicted to come from stromal, 23.2% immune, 16.8% neuronal, and 23.2% epithelial). Recombinant human proteins (Supplementary Data 13) were applied to differentiated cytomix-treated colon organoids from 3 donors in an arrayed format followed by multiplexed scRNA-seq to assess transcriptional responses (Fig. 3A).

To map the overall landscape of relationships between perturbations, we calculated the cosine similarity between all perturbations (Fig. 3B). Ligands in clusters 3–8 (Low Response) display similarities to cytomix-only controls, suggesting that the majority of these perturbations did not elicit an identifiable response over the injury itself. Perturbations from the same family generally clustered together with high similarity at the pseudo-bulked level (EGF-family ligands, TGFβ-superfamily ligands, pro-inflammatory molecules IFNε, TNFɑ, IFNɣ), demonstrating the ability to recover known ligand pathway redundancies. We also uncovered divergences in activity across ligand families, including NRG1 clustering with SLIT3 and NT-3 instead of other EGF-family ligands, potentially due to EGF-like domains within SLIT3, and TRAIL (TNFSF10)/TWEAK (TNFSF12) clustering with the TGFβ-superfamily rather than RANKL (TNFSF11) and TNFɑ, suggesting disparate responses in the context of the inflamed epithelium. Interestingly, uninjured (no cytomix) organoids clustered with members of the TGFβ-superfamily, including BMP6 and GFD11, suggesting that BMP treatment on inflamed organoids may elicit a transcriptional response that drives higher similarity to non-inflamed organoids. Highly correlated perturbations are also preserved at the cell-type-specific level indicating similarities are largely not driven by cell-type composition differences.

To determine the genes driving these relationships, we performed DEG analysis for every perturbation vs. cytomix-only controls in a cell-type-specific manner (Supplementary Fig. 7A–D, Supplementary Data 7). 25/80 perturbations elicited a response with more than 1 DEG (false discovery rate (FDR) < 0.05, Fig. 3C). In particular, members of the BMP, TGFβ, TNF, and IL families exhibited strong transcriptional responses with high numbers of DEGs (>25 based on p-adj < 0.05, abs(log2FC) > 1) and donor reproducibility (Fig. 3C). We recovered known biology as ligands upregulated canonical target genes, including BMP2/6 (ID1, ID3), TGFβ1/2/3 (SMAD3, CDKN1A), and IL22 (DMBT1, MUC1)51. We additionally uncovered less well-known responses which may be specific to an inflamed context, including BMP2/6 upregulating genes associated with GP IBD2 (CKB, SELENBP1, C15orf48) suggesting a role in maintaining metabolic homeostasis, TGFβ1/2/3 upregulating interferon response genes (SP100, CX3CL1) implicating their roles in epithelial response to inflammation, and TRAIL downregulating pro-apoptotic genes (BNIP3L, FAM162A, BNIP3) suggesting an anti-inflammatory role (Supplementary Fig. 7E).

To gain a global view of shared regulatory modules across perturbation clusters, we organized DEGs into gene programs and determined the average effect size of each GP for each perturbation versus cytomix-only controls (Fig. 3D, Supplementary Fig. 6C, Supplementary Data 8). GPs with the highest effect size continued to be centered around the TGFβ-superfamily and Inflammatory perturbation clusters, confirming the transcriptional impact of these perturbations. Not only did we observe unexpected similarities in GP activity between dissimilar ligand families, but also divergent GP activation between similar ones. For instance, BMP2/4/6 and TGFβ1/2 downregulated GPs 1 and 5, while BTC, HBEGF, and EGF upregulated these GPs, indicating a shared ability of both ligand families to modulate these programs. Conversely, GP6 was upregulated by BMP2/6 but downregulated by BMP4.

We next leveraged the single cell resolution of our atlas to investigate cell-type-specific response across secreted ligands. Some perturbations led to the emergence of new cell populations, such as the M cell-like population in RANKL-treated organoids (Fig. 4A)54. Other perturbations elicited drastic changes in cell-type composition. BMP2, 4, and 6 led to a decrease in the proportion of stem cells and increase in colonocytes, consistent with BMP’s role in restricting stem cell hyperproliferation13,55 (Fig. 4B, Supplementary Fig. 8). On the other hand, inflammatory factors like IFNɣ and TNFɑ led to a decrease in colonocyte proportion (Fig. 4B, C), in line with observations in human diseased vs. healthy IBD patient samples4. EGF, HBEGF, IFNɑ, RANKL, and overall cytomix stimulation decreased goblet cell proportions while IL22, BMP6, and TRAIL (TNFSF10) caused an increase. Not only has recent evidence suggested that increased goblet cell abundance is associated with clinical remission following treatment, goblet cells in TRAIL-treated organoids displayed increased expression of remission-associated goblet markers (MUC2, SPINK4, LCN2) (Supplementary Data 7)6.

To further understand how DEGs are co-regulated in a cell-type specific manner, we calculated effect sizes of previously determined GPs (Fig. 3D) across individual cell types for select high impact perturbations (Fig. 4D, Supplementary Data 11). Cytomix treatment itself exerted the strongest effects, as no-cytomix organoids downregulated GP20 (immune response to cytokines: IFITM3, IFI6, DDX60) and GP13 (Neutrophil chemotaxis, response to cytokine: CXCL9, CXCL10, CXCL11), and upregulated GP2 (ribosome biogenesis: RPS12, RPS28, RPL39) across all cell types. GP7, related to organ growth and cell migration (EDN1, ANXA3, WWC2) was downregulated most strongly in differentiated cells with BMP2, 6, and 4 but not other perturbations. Interestingly, besides GP7, additional responses to BMP4 appeared to be more transcriptionally distinct from BMP2/6, suggesting divergences in the BMP family ligands (Supplementary Fig. 7A–D). GP6, related to proliferation and regulation of wound healing (PTPRK, SH3RF2, HEPH) exhibited opposing effects in BMP2/6 compared to EGF-family ligands, particularly in stem/TA cells. Notably, treatment with members of the TGFβ-superfamily, particularly BMPs, resulted in similar trends to no-cytomix organoids, suggesting that specific aspects of the non-inflamed epithelium can be recapitulated by administration of TGFβ-superfamily ligands in the presence of cytomix. Specifically, GP1 (Ion transport and cell adhesion: PLCB4, ITPR2, SKAP1), is downregulated in stem and TA cells in both no-cytomix and BMP2/6/TGFβ ligands and upregulated in TA cells for pro-inflammatory perturbations TNFɑ and IFNɣ. In contrast, GP12 (epithelial structure maintenance and downregulation of anoikis: FABP1, LGALS3, S100A14) is upregulated in stem and/or TA cells in both no-cytomix and BMP2/6/TGFβ ligands and downregulated in IL13, IL4, and pro-inflammatory perturbations (Fig. 4D, E).

To probe the functional relevance of GPs elicited in organoids, we searched for their expression in human IBD tissue4,25 (Fig. 5A). We focused on GP1, 12, and 20 as they show strong effect sizes in cytomix treatment and TNFɑ/IFNɣ individually, suggesting potential disease-relevance. These GPs additionally display large effect sizes in organoid upon stimulation with perturbations in the TGFβ-superfamily and IL4/13 (GP1-downregulated by TGFβ-superfamily; GP20 - downregulated by BMP6 and IL4; GP12 - upregulated by TGFβ-superfamily, downregulated by cytomix). Notably, many genes in GP1, 12, and GP20 showed a decrease in healthy vs. IBD epithelium across two scRNA-seq datasets (Fig. 5A, Supplementary Data 10). To further examine the localization of these GPs in human tissue, we performed spatial transcriptomics (CosMx) on healthy (n=4) vs UC (n=16) patient samples (Supplementary Data 9) (Fig. 5B). GP1 was excluded from this analysis due to low panel coverage. Not only were GP12 and GP20 epithelial-enriched, they were also upregulated in UC vs. normal epithelium (Fig. 5B, C). Altogether, this suggests that organoid-derived GPs are capable of modeling physiologically relevant responses.

Due to the cyclic nature of damage and repair in IBD, it can be difficult to interpret whether dysregulated gene signatures in human tissues are related to pathogenesis vs. regeneration. However, leveraging organoid-derived GPs, particularly from perturbations with known biology (e.g. TNFɑ and IFNɣ are known to drive inflammation56,57), we can shed light on potential directionality. For instance, based on the expression patterns of GP1, 12, and 20 in organoids treated with inflammatory cytokines, we infer that upregulation of GP1 and 20 (up with cytomix) in human tissues may play pro-inflammatory roles, while upregulation of GP12 (down with cytomix) may be protective. Based on expression patterns of such GPs, we can further infer that members of the TGFβ-superfamily may play protective roles in the inflamed epithelium58–62.

Finally, we checked the expression of IBD-related GPs IBD1 and IBD2 (Fig. 1) across perturbations (Fig. 5D, E, Supplementary Data 11). Inflammation-associated IBD1 was slightly downregulated in TA cells from BMP6, IL13, and IL4-treated organoids60. We observed a stronger response for metabolism-associated IBD2, which was upregulated across the majority of cell types in no-cytomix along with BMP2/6, and in stem cells for TGFβ1/2/3 (Fig. 5D, E). This further suggests that some members of the TGFβ-superfamily family of ligands as well as IL4/13 may recapitulate aspects of transcriptional states found in the uninflamed epithelium.

Methods
See Supplementary Data 12 for more information on reagents used in this study.

Ethics statement
All human tissue samples were supplied by Donor Network West, who received IRB approval of research, appropriate informed consent of all subjects contributing biological materials, and all other authorizations, consents, or permissions as necessary for the transfer and use of the biological materials for research at Genentech. Patients/human donors were not compensated for their donation of tissues.

Generation of primary human colon organoids
Sections of human transverse colon were provided by Donor Network West from deceased patients with no history of intestinal disease. For the perturbation atlas, organoids from 2 female donors and 1 male donor were used. To clean the tissue and isolate fresh crypts for organoid generation, colon tissue was first sliced open lengthwise and wiped 3–4 times with gauze (Fisher Scientific) to remove waste and mucus. The tissue was then cut into small sections of ~2”×2” and placed into a 15 cm dish filled with PBS++ (PBS supplemented with 10 mM HEPES (Thermo Fisher Scientific), 1X GlutaMAX (Thermo Fisher Scientific), 1X Pen-Strep (Thermo Fisher Scientific), and 2.5 μg/mL Amphotericin B (Thermo Fisher Scientific)). The epithelium was then minced from tissue segments using scalpels (Fisher Scientific) and collected into 50 mL conical tubes. Minced epithelium was washed 3 times with DMEM++ (Advanced DMEM/F12 (Thermo Fisher Scientific) supplemented with 10 mM HEPES, 1X GlutaMAX, 1X Pen-Strep, and 2.5 μg/mL Amphotericin B), and centrifuged at 4°C, 500×g for 3 min between each wash. After the final wash, tissue was resuspended in 25 mL of 2.5 mM EDTA in PBS and placed in an end-over-end rotator at 37°C until crypts were visible, around 12 min. Following digestion, the EDTA digestion solution was topped up with 20 mL of DMEM++, filtered through sterile gauze, and centrifuged for 3 min at 4°C, 500×g. The pellet was then resuspended in 10 mL PBS++, filtered through a 100 μm cell strainer, and centrifuged for 3 min at 4°C, 500. The tissue was then ready for plating in Matrigel (Corning) by resuspending in cold Matrigel and plating into 30 μL Matrigel droplets, at a density of 30–50 crypts per droplet. Once Matrigel had solidified for 30 min at 37°C, crypts were fed with organoid growth media65: Advanced DMEM/F12 supplemented with 1X Pen-Strep, 10 mM HEPES, 1X GlutaMAX, 1X B-27 Supplement (Thermo Fisher Scientific), 0.2 nM Wnt Surrogate protein (ImmunoPrecise Antibodies), 5 nM Gastrin-1 (Tocris), 1.25 mM N-Acetyl-L-Cysteine (Sigma-Aldrich), 10 mM Nicotinamide (Sigma-Aldrich), 100 μg/mL Primocin (InvivoGen), 500 nM A83-01 (Tocris), 50 ng/mL recombinant human EGF (R&D Systems), 100 ng/mL recombinant human Noggin (R&D Systems), 250 ng/mL recombinant human R-spondin 3 (Roche), 100 ng/mL recombinant human IGF1 (BioLegend), 50 ng/mL recombinant human FGF2 (PeproTech).

Organoid maintenance and passaging
After the initial organoid plating, cultures were maintained in Cultrex Basement Membrane ECM (70% Cultrex-PBS) (R&D Systems) domes. Organoids were cultured every 5–7 days as needed based on size and morphology. To passage, organoids were mechanically dislodged from Cultrex domes using a P1000 pipette tip and resuspended in cold 10 mM EDTA with 10 μM Y-27632 (STEMCELL Technologies) in a 15 mL conical tube. The tube containing organoids was then placed on an end-over-end rotator at 4°C for 30 min and then centrifuged for 5 min at 500×g, 4°C. Organoids were then resuspended in 1X Trypsin (Fisher Scientific) and placed in an end-over-end rotator for 6 min at 37°C. Following dissociation, Trypsin was neutralized with DMEM containing 10% FBS, and organoids were centrifuged for 5 min at 500×g, 4°C. Media was aspirated and organoids were washed once with PBS, followed by centrifugation for 5 min at 300×g, 4°C. PBS was aspirated and organoids were resuspended in 70% Cultrex with PBS, and plated in 30 uL droplets in 6 well plates. The plates were placed upside down in the incubator (37°C, 5% CO2) for 30 min, after which organoid growth media with 10 μM Y-27632 was added. Media was changed every other day to organoid growth media.

Organoid culture in microwell plates
Experiments involving microwell plates utilized 500 μm Gri3D plates (SUN bioscience). Prior to organoid seeding, storage buffer was aspirated from each well and 150 μL of Advanced DMEM/F12 was added through the pipetting port. The plate was placed in the tissue culture incubator for at least 30 min. Organoids that had been maintained in Cultrex domes for at least 5 days after passaging were removed from the ECM following the passaging protocol as above, but dissociated for 12 min in Trypsin rather than 6, or until organoids were largely single cells. Trypsin was neutralized and organoids were washed as in the passaging protocol. Following the PBS wash, organoids were resuspended in organoid growth media without Nicotinamide: Advanced DMEM/F12 supplemented with 1X Pen-Strep, 10 mM HEPES, 1X GlutaMAX, 1X B-27 Supplement, 0.2 nM Wnt Surrogate protein, 5 nM Gastrin-1, 1.25 mM N-Acetyl-L-Cysteine, 100 μg/mL Primocin, 500 nM A 83-01, 50 ng/mL recombinant human EGF, 100 ng/mL recombinant human Noggin, 250 ng/mL recombinant human R-spondin 3, 100 ng/mL recombinant human IGF1, 50 ng/mL recombinant human FGF2, and with 10 μM Y-27632. Dissociated organoids were then counted using a hemocytometer and seeded into the microwell plates at a density of 200 cells per microwell, in 50 μL of media per well. The plates were placed in the incubator for 30–60 min to allow the cells to settle, after which another 100 μL of media containing 2% Matrigel was added per well. It is critical that Matrigel is added to the media while cold. On the day after seeding, media was removed and changed to growth media without Y-27632. On day 3 after seeding, media was removed and organoids were washed twice with Advanced DMEM/F12. After washing, media was replaced with organoid media without EGF to promote differentiation. Media was changed every 1–2 days, and cultures were maintained in the plates for up to 10 days before collection.

Injury models in organoids
For experiments involving irradiation, DSS, and cytomix, organoids were cultured as above in microwell plates for 5 days. On the 5th day after seeding, injury conditions were applied. For irradiation studies, the plates were placed on the center of the rotating table of the irradiator (J.L. Shepherd and Associates, Mark I series, model 68A, Cesium (137Cs) gamma irradiator) for the time periods equivalent to either 6Gy or 12Gy, and then the organoids were further cultured in the incubator. For DSS treatment, DSS (MP Biomedicals) was dissolved in organoid media (without EGF or Nicotinamide), and added to organoid wells overnight for 16 h, after which organoids were washed twice with Advanced DMEM/F12 and then replaced with organoid media. The low and high doses considered were 0.1% DSS and 0.5% (wt/vol) DSS in media. For cytomix treatment, human recombinant TNFɑ (R&D Systems), IFNɣ (R&D systems), and IL-1β (PeproTech) were dissolved in organoid media (without EGF or Nicotinamide) for 8 hours, after which organoids were washed twice with Advanced DMEM/F12 and then replaced with organoid media. The low and high doses considered were 10 ng/mL TNFɑ, 10 ng/mL IFNɣ, 1 ng/mL IL1β and 50 ng/mL TNFɑ, 50 ng/mL IFNɣ, 5 ng/mL IL1β. The high and low injury conditions were selected on pilot data to evaluate dose-response (Supplementary Fig. 2A). Experiments were conducted on the following day after injury (scRNA-seq, imaging) to assess response to injury.

FITC-dextran and apoptosis measurements in organoids
Organoid permeability and apoptosis measurements were conducted on one day post-injury, and FITC-dextran was used to assess barrier integrity34. Briefly, organoids were incubated in a solution of 4kDA FITC-Dextran (MedChemExpress) diluted to 2 mg/mL in PBS for 5 min in μ-slides (Ibidi) to allow organoids to settle to the bottom of the slide chamber. Organoids were then imaged on a Nikon spinning disc confocal to collect fluorescence and DIC images, and the percentage of intact organoids was calculated as the number of organoids per condition/donor that excluded FITC-dextran from the lumen divided by the total number of organoids. For apoptosis quantification, organoids were incubated in CellEvent Caspase-3/7 Detection Reagent (Thermo Fisher) diluted in media according to manufacturer protocols for 30 min in the microwell plates. Organoids were then washed with PBS and imaged on a Nikon spinning disc confocal to collect fluorescence and DIC images. The % area of caspase was quantified using FIJI66 as the area of organoid stained positive for Caspase-3/7 divided by the total area of the organoid from DIC images. Sample sizes for permeability measurements in Fig. 1C were taken from n=3 independent batches of organoids, each containing at least 16 organoids per condition (exact numbers shown in Source Data). Sample sizes for Caspase-3/7 measurements in Fig. 1C are control=12, DSS_low=11, DSs_high=10, cytomix_low=13, cytomix_high=10, irradiation_low=10, and irradiation_high=10, where numbers represent distinct organoids.

Organoid whole-mount immunofluorescence
For immunofluorescent images, organoids were removed from the microwell plates and fixed in 4% paraformaldehyde (Thermo Fisher) overnight at 4°C in 1.5 mL Eppendorf tubes on a rocker. On the next day, organoids were washed three times with PBS and resuspended in permeabilization/block buffer (5% normal donkey serum (Abcam), 0.5% Triton X-100 (Thermo Fisher) in PBS)67 overnight at 4°C on a rocker. Organoids were then resuspended in primary antibody diluted in organoid wash buffer (0.1% Triton X-100, 0.2% BSA (Miltenyi)) overnight at 4°C on a rocker. On the following day, organoids were washed three times (2 h per wash at 4°C on a rocker) in organoid wash buffer, and then resuspended in secondary antibody with DAPI diluted in organoid wash buffer. On the next day, organoids were washed three times with PBS (2 h per wash at 4°C on a rocker), and then placed in a 96-well glass bottom plate. PBS was carefully removed and organoids were resuspended in 60 μL of CUBIC-R (TCI). After allowing organoids to settle, imaging was conducted on a Leica STED. Brightness and contrast adjustments were performed using Fiji66, and adjustments were made uniformly across images. Primary antibodies used are Mouse monoclonal anti-E-cadherin 1:500 (BD Transduction Laboratories, Cat# 610181), Goat polyclonal anti-E-cadherin 1:500 (R&D Systems, Cat# AF748), Rabbit monoclonal anti-KI67 1:200 (Thermo Fisher, Cat# RM-9106-S1), Mouse monoclonal anti-MUC2 1:200 (Agilent DAKO, Cat# IR65861-2), and Rabit polyclonal anti-CHGA 1:200 (Novus Biologicals, Cat# NB120-15160). Secondary antibodies are purchased from Jackson ImmunoResearch and used at a dilution of 1:500: Donkey anti-mouse 488 (Cat# 715-545-150), Donkey anti-mouse Cy3 (Cat# 715-165-150), Donkey anti-mouse 647, Cat# 715-605-150), Donkey anti-rabbit 488 (Cat# 711-545-152), Donkey anti-rabbit Cy3 (Cat# 711-165-152), Donkey anti-rabbit 647 (Cat# 711-605-152), Donkey anti-goat 488 (Cat# 705-545-147), Donkey anti-goat Cy3 (Cat# 705-165-147), and Donkey anti-goat 647 (Cat# 705-605-147).

Single cell preparation of organoids for single cell RNA sequencing
For scRNA-seq experiments, organoids were dissociated using the Neural Dissociation-P kit (Miltenyi). Organoids were pipetted out of their microwell plates and into a deep-well 96 well plate (Thermo Fisher) for dissociation. Enzyme mix 1 was prepared according to manufacturer protocols with the addition of Y-27632 for a final volume of 1.5 mL per well. Organoids in mix 1 were then placed at 10°C for 10 min, after which enzyme mix 2 was added. Organoids were kept at 10°C throughout dissociation for up to 60 min, and were removed every 10 min to be pipette mixed with a P1000 pipette to aid in dissociation. Once dissociation was complete, organoids were centrifuged for 3 min, 500×g, 4°C and washed twice with PBS. For pooled experiments involving hashing (perturbation experiment, pooled injury experiments), organoids were resuspended in 100 μL cell staining buffer (2% BSA, 0.01% Tween in PBS) and then blocked with 5 uL human TruStain FCX blocking buffer (BioLegend) for 10 min on ice. 1 μL of a unique cell hashing antibody (TotalSeqA: BioLegend) was added per condition for 20 min on ice, and cells were pipette mixed after 10 min. Organoids were then washed three times with cell staining buffer on ice. After the final wash, organoids were resuspended in 0.04% BSA in PBS, filtered through a 70 μm Flowmi cell strainer (Sigma Aldrich) and loaded onto a 10X Chromium X according to manufacturer instructions. For organoid and tissue characterization and pooled injury experiments, the standard 10x 3’ v3.1 kit was used. For some injury experiments, the 10x 3’ low throughput kit was used to compare multiple injuries from the same experiment. For the perturbation dataset, the 10x 3’ high throughput kit was used. Gene expression libraries were prepared according to 10X protocols, and HTO libraries were prepared according to published protocols68.

scRNA-seq data pre-processing
Libraries were sequenced on a NextSeq 2000 or NovaSeq 6000 (for HT libraries) targeting 30,000 reads/cell for GEX libraries and 2000 reads/cell for HTO libraries. Sequenced reads were aligned using the human reference genome GRCh38-2020-A using the Cell Ranger pipeline (10x Genomics, v7.1.0). scRNA-seq data was processed using CellRanger and resulting count matrices were analyzed using (v5.0.2) in Rstudio (v4.2.1). For experiments involving hashing, cells were demultiplexed using the HTO Demux function in R. For experiments involving genetic demultiplexing as in the perturbation dataset, the souporcell algorithm69 was used to distinguish between organoids from unique donors. Doublets were removed from data using the DoubletFinder algorithm (v2.0.4)70. When combining experiments from multiple donors or timepoints, data was batch-corrected using the Harmony algorithm (v0.1.1)71.

Secreted niche factor perturbation screen
Human colon organoids from 3 independent donors were expanded in Cultrex domes and then seeded in microwell plates as described in the above section. On the 5th day post-seeding, low-dose cytomix (10 ng/mL TNFɑ, 10 ng/mL IFNɣ, 1 ng/mL IL1β) was applied to the organoids in a media change with 0.5% Matrigel. At the time of cytomix addition, human recombinant proteins corresponding to the 79 selected perturbations (Supplementary Data 13) were applied to each well in an arrayed format. Cytokine concentrations were selected to be within or above the EC50 listed from the vendor product information in order to increase the chances of observing a transcriptional effect, which for most compounds was 100 ng/mL. For our pro-inflammatory cytokines TNFɑ and IFNɣ, we have observed evidence of toxicity and cell death in preliminary experiments at 100 ng/mL, so we reduced the concentrations of these molecules down to 50 ng/mL. Additionally, we increased the concentrations for some factors (ANGPTL1, ANGPTL4, AREG, BDNF, BMP8B, CT-1, CTSG, F2, FGF1, FGF18, FGF7, FGF9, NT-4, POSTN, PPY, TWEAK, VIP) up to 250 ng/mL based on preliminary data in our lab showing that these compounds did not elicit strong transcriptional effects at 100 ng/mL. After 8 h, cytomix was removed from the organoids by washing twice with Advanced DMEM/F12 for 30 min in the incubator. After washing, fresh organoid media with 0.5% Matrigel and human recombinant proteins was added to each well (Supplementary Data 13). Media was changed every 1–2 days with fresh organoid media containing human recombinant proteins, and organoids were collected for scRNA-seq on day 8 post-seeding. Dissociated organoids were loaded onto 4 3’ HT channels (10X), and for each independent channel one no-cytomix control as well as 3 cytomix-only controls were included.

Spatial transcriptomics (CosMx) data collection for human colon tissue
Formalin-fixed, paraffin-embedded (FFPE) tissue blocks were sectioned at 5 μm thickness using a microtome and processed according to Nanostring’s protocol for “CosMx SMI Semi-Automated Slide Preparation for RNA Assays”. Briefly, the slides were placed in an oven at 60°C overnight to facilitate paraffin removal and sample conditioning. Deparaffinization and antigen retrieval were performed using the Leica Bond RX system. Slides were incubated with ER1 (AR9961, Leica) at 100°C for 20 min to retrieve antigens. Subsequently, proteinase K (AM2546, Thermo Fisher) was applied at a concentration of 3 μg/mL for 15 min to enhance tissue permeability and facilitate probe hybridization. Fiducial markers were applied for 5 min to aid in downstream imaging. Afterwards, slides were incubated with 100 mM NHS-acetate for 15 min to block nonspecific binding sites, followed by washing with 2X SSC buffer. Post-fixation was carried out using 10% neutral buffered formalin (16004-122, VWR) for 5 min. Subsequently, slides were incubated with 100 mM NHS-acetate for 15 min to block nonspecific binding sites, followed by washing with 2X SSC buffer. Denatured CosMx RNA probe mixes were applied to the slides. Hybridization was performed overnight under nuclease-free conditions (RNase away and DEPC treated water) in a hybridization oven (HybEZ, ACD) at 37°C. After hybridization, stringent washes were performed twice using a stringent solution (50% formamide in 2X SSC buffer solution) at 37°C for 25 min each. Unbound probes were removed and background signal reduced by washes with 2X SSC buffer solution. Slides were incubated with nuclear stain for 15 min, followed by application of cell segmentation mix (CD298/B2M) with supplemental markers (PanCK/CD45) for 1 h. Slides were washed with 1X PBS to remove excess staining reagents and flow cells were assembled to the slides. To prevent tissue desiccation, 2X SSC buffer was pipetted into the flow cell ports.

The prepared slides and the CosMx RNA imaging 1000-plex tray were loaded onto the CosMx Spatial Molecular Imaging (SMI) platform. The flow cell setting was created for pre-bleaching profile using Colorectal Tissue (configuration A) and the cell segmentation profile used was Human (configuration A). Scan area and fields of view (FOV) were selected for each slide and data acquired during the cycling phase on the CosMx (SMI) platform, capturing high-resolution, spatially resolved RNA expression profiles. During the CosMx cycling, data was automatically transferred to the AtoMx platform for further analysis.

Cell type annotation in organoids
Organoid scRNA-seq data was quality controlled, filtered and annotated using Seurat (v5.0.1)72. Cell-type annotation was performed by applying FindTransferAnchors and TransferData using our matched colonic tissue reference atlas (Supplementary Fig. 1), and then manually cleaning up annotations based on expression of canonical cell type markers. The human colon organoid tissue atlas was first annotated by integrating our colon data with published atlases containing healthy colon tissue4,25,73 using scANVI to harmonize cell type labels in a consistent manner, thus generating an integrated atlas74.

Derivation of colon epithelial IBD signatures from published atlases
scRNA-seq data downloaded from public datasets (Smillie et al., Kong et al., and Parikh et al.)4,25,48 were imported and cells with greater than 30% mitochondrial gene content and fewer than 200 genes expressed were removed from all datasets for downstream analysis. Cell-type annotation was re-performed by projecting each query dataset onto our human colon integrated atlas (using Seurat functions FindTransferAnchors and TransferData). Filtered RNA counts matrices were converted to h5ad format via the scanpy (v1.9.3) for input into consensus non-negative matrix factorization (cNMF)75. The cNMF package (v1.3.4) was utilized to identify gene modules with a number of components (k) set to 10, a local density threshold of 0.2, and a maximum of 10 iterations. Metamodules were identified by calculating the Spearman’s correlation using the gene spectra score space of cNMF modules derived independently from Smillie and Kong datasets. The 50 genes with greatest spectra scores from the modules within each meta-module were identified and the unions of these gene sets were the IBD1 and IBD2 colon epithelial IBD signatures.

For both the effect-size dotplots and the boxplots visualizing enrichment of metamodules in patients with inflamed vs non-inflamed conditions, IBD1/IBD2 signature scores were computed using the AddModuleScore function in Seurat. The means of module scores for each patient were then calculated, with n=20 Kong et al. IBD, n=16 Kong et al. healthy, n=6 Parikh et al. IBD, n=3 Parikh et al. healthy, n=18 Smillie et al. IBD, and n=12 Kong et al. healthy. Effect sizes were evaluated using Cohen’s d formula for inflamed and non-inflamed groups for each cell type.

The predictive power of IBD1 and IBD2 for disease status was assessed by training a generalized linear model (GLM) on the Kong et al. dataset25, with 70% of the data used for training and testing on the remaining 30%, and finally validating in Smillie et al. dataset4.

Cosine similarity, DEG, gene program and effect size analyses
Cosine similarities between perturbations was calculated on the Harmony donor-corrected PCA space (first 50 PCs). An average PC coordinate was calculated for each perturbation. To calculate the donor-donor correlation, average PC coordinates were calculated on a donor basis and cosine similarities were calculated for each donor-perturbation pairwise comparison. Significance was met if the average donor-donor correlation for a single perturbation was >95% percentile of a permuted distribution (rho=0.406).

For the injury dataset, DEGs were calculated for each injury vs. control for each cell type, by fitting the data to a gamma-possion model via the glmGamPoi package using only the top 5000 highly variable genes and a randomly downsampled 300 cells per cell type per injury76. Log FC values were shrunken. For the perturbation dataset, DEGs were calculated by pseudobulking counts for each gene on a donor basis and by cell types. DEseq2 was used to compare each perturbation to the untreated control, with the design testing for perturbation while controlling for donor77. DEGs were designated as padj < 0.05 and log2FC > 0.2. The average number of DEGs for each perturbation was calculated as the average across each major cell type for each perturbation.

To group perturbation-elicited DEGs into biologically interpretable gene programs, we ran cNMF78 on all DEGs in each dataset with a number of components (k) set to 20, a local density threshold of 0.2, and a maximum of 10 iterations. The resulting module usages were utilized for analyses comparing GP expression. Activity vs. identity modules were manually selected based on expression of canonical cell type-specific markers and restriction of module expression to particular cell types, and only cell activity modules were considered for further evaluation.

In the 83-perturbation data set, the effect size (Cohen’s d) was calculated using the module usage scores (from cNMF), and each perturbation (3 donors, 1 well per donor) was only compared against the controls (3 donors, 3 wells per donor) that were from the same cDNA/library pool, to control for pool effects. The standard deviation and mean was taken across 3 donors for each perturbation. The p-value for each effect size was calculated using a two-sided t-test on the average module usage scores on a well/donor-basis, and FDR corrected.

Effect sizes of IBD1 and IBD2 were calculated by approximating the average expression of genes in each GP with AddModuleScore in Seurat, considering the top 50 genes in each GP based on gene score. For the injury dataset, the mean/standard deviation was taken across metacells computed across the dataset using the hdWGCNA package79 in R, with k=20 and max_shared=4.

Gene Ontology analysis: pathway enrichment of selected gene programs
GO analysis was performed using the enrichR package (v3.1)80 in R, specifically querying the GO biological processes database from 2023. The top 50 genes from each GP based on gene spectra score were selected for inclusion in GO analysis.

Selection of secreted factors for perturbation screen through ligand-receptor analysis
The perturbation space was established by mining through published scRNA-seq atlases of human colon from development (Yu et al., Elmentaite et al.)26,27, adult homeostasis (Elmentaite et al., Kong et al.)25,26, and IBD (Kong et al.)25. Data from each atlas was considered individually and filtered to only include samples from the colon. Ligand-receptor interactions were assessed in each atlas using CellChat (v1.6.1)81, with trim=0.02 and min.cells=10. From each dataset, perturbations signaling to the epithelium were selected by setting targets.use to each epithelial cell type as annotated in the original dataset. After running CellChat on each dataset, perturbation lists were combined into a spreadsheet documenting each predicted signaling interaction and the probability of this interaction (Supplementary Data 5), and prediction probabilities were averaged across all predicted interactions across datasets. The perturbation space was cut down to 79 to fit onto one 96-well plate (arrayed format) including controls by removing low-probability perturbations, perturbations exclusively signaling from immune cells, and preferentially selecting or adding in ligands from the same family or with known biological function (IFNɑ, IL1ɑ, BMP3, DKK1, EdIL3, IFNE, IL4, SLIT3, ISLR, IL22).

Calculation of log FC and p values in external datasets
For external IBD datasets4,25, DEGs were calculated by pseudobulking counts for each gene on a donor basis and by cell types. DESeq2 was used to compare healthy vs. diseased cells, with the design testing for disease state while controlling for donor. DESeq2 models the data using a negative binomial distribution, shrinkage of log FC values was applied using the apeglm method, and significance was calculated using a two-sided Wald test.

Cell type composition analysis
Changes in cell type composition upon perturbations compared to the Control were determined using a beta-binomial regression framework via the DCATS package82, with n=3 donors per perturbation. Results were additionally confirmed via differential abundance testing via MILO83 by considering donor and pool as covariates.

Spatial transcriptomics (CosMx) analysis
Analysis of spatial transcriptomics data generated by the CosMx SMI 1000-plex protocol was performed in two python workflows following standard quality control procedures output by the AtoMx segmentation and transcript assignment: 1) ensemble variational inference and label transfer and 2) gene set analysis and visualization.

First, ensemble variational inference and label transfer was used to annotate segmented single-cell spatial transcriptomes at a high resolution using a publicly available scRNA-seq dataset as a reference39. We performed 10 replicate runs of scVI followed by scANVI as a Papermill workflow on an Nvidia Tesla V100 32GB GPU. scVI parameters were as follows: n_latent=32; max_epochs_scvi=100; dropoutrate=0.3 (for resection samples) and 0.1 (for TMA samples); and batch_key=“Modality”, where we encoded “Reference” or “Spatial query” corresponding to the appropriate datasets. scANVI parameters were as follows: n_samples_per_label: 1000. In order to ensemble the results of all replicate runs we applied a classification entropy filter followed by consensus voting. Each replicate has an associated cell ID by cell type classification probability matrix, where the cell type classification entropy was calculated cell-wise and thresholded by finding the lowest local minimum by examining the histogram. Cells with classification entropies above this threshold are considered “Unknowns”, interpreted as mixed single-cell transcriptomes that could not be confidently assigned to a single cell type. Next, consensus voting is performed and the final annotation is assigned by majority voting given a proportional threshold of 0.5, where classifications below this number are again assigned as unknowns. For cell type-oriented analyses these unknowns were excluded. Further lineage-specific cell type refinements were performed through high-resolution leiden clustering on the output scANVI embedding.

Second, spatial gene set analysis was performed using the scanpy and squidpy python libraries. These data were normalized by library size and log1p transformed before scoring using the scanpy score_genes function with default parameters on each single cell. Gene set scores were calculated only on the intersection of the 1000-plex CosMx panel. For GP analyses in the context of the TMA samples, this scoring procedure was performed and averaged per biological replicate and per cell type. Mann–Whitney U tests were performed to assess statistical significance. For GP analyses in the context of the resected normal and UC samples, the same scoring procedure was performed; to visualize, the resulting scores were binned into 20 quantiles and plotted by their assigned 2D spatial coordinates.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This study presents a valuable single-cell and spatial transcriptomic atlas of human colonic epithelial responses to inflammatory stimuli and niche factor perturbations. The manuscript offers a comprehensive resource that maps epithelial responses to microenvironmental cues in an IBD-relevant organoid model, effectively linking in vitro findings to human disease through scRNA-seq and spatial analyses. Non-biased and comprehensive study is entirely valuable data source for proceeding the related field.\nOne comment is the connection with fetal reversion. or about recent. The authors identify a CXCL8⁺AREG⁺ANXA1⁺ population in dissociation-induced organoids (Fig. 1F), which resembles previously described revival stem cells (revSCs) as reported in several seminal studies (Yui et al., Cell Stem Cell, 2018; Nusse et al., Cell Stem Cell, 2018; Ayyaz et al., Nature, 2019; Klementaitė et al., Nature, 2020). However, the potential presence and physiological relevance of such revSC-like populations in IBD patient tissues are not explored. Given that the authors have access to scRNA-seq and CosMx spatial transcriptomic data from ulcerative colitis patients, they are in a strong position to assess the expression and localization of canonical revSC markers such as ANXA1, TACSTD2, and CLU in vivo. Incorporating such an analysis would significantly strengthen the interpretation of the regenerative programs described in the organoid system.\nIn addition, the gene programs (GPs) extracted in this study are mainly based on in vitro data generated under a specific matrix condition: all organoids were cultured in 2% Matrigel suspension across injury and perturbation conditions. As shown in Figure 1B, the epithelial cells are basically polarized as written in the text. While this uniform setup facilitates comparative analysis, it is important to recognize that matrix composition and mechanical context can profoundly influence epithelial cell states, particularly under regenerative or inflammatory conditions. For example, Yui et al. (2018) showed that revSC signatures emerge preferentially in collagen-based cultures with the features of depolarization of epithelial cells. Moreover, the authors themselves highlight “cell–matrix adhesion” as a defining component of GP6 (Fig. 4D), suggesting ECM-dependent transcriptional responses.\nI do not think that additional experiments are required. However, I would recommend that the authors consider the following:\n1. Acknowledge the limitations of their matrix system in shaping regenerative cell states. For example, human colon collagen-based organoids have been shown to be enriched in several CXCL genes (Kobayashi S. et al., Inflammation and Regeneration, 2022), which may not be fully captured in Matrigel-based systems.\n2. Discuss the possibility that certain regenerative programs—such as revSC-like states—may be underrepresented or modulated by the current culture context. Alternatively, if such states are embedded within specific GPs identified in this study, it would be valuable to relate those programs to the known revSC phenotype and situate the findings within the broader context of intestinal inflammation and epithelial regeneration. Even though the perturbation screen is non-biased, the resulting transcriptional programs show substantial overlap with features of previously described revSC states.\n3. Include key references on revSC biology (Yui et al., 2018; Nusse et al., 2018; Ayyaz et al., 2019; Klementaitė et al., 2020 etc., if any others important) to better contextualize the regenerative populations described here and enhance the biological relevance of the findings.\nIncorporating these points would improve mainly the Discussion section by providing a more comprehensive and physiologically informed interpretation of the epithelial regenerative responses uncovered in this study.", "reviewer_2": "Capeling MM et al. present single-cell transcriptomic data of human colonic organoids subjected to various cytokine stimulations. They also integrate their single-cell dataset with spatial transcriptomic data from inflammatory bowel disease (IBD) tissues. This dataset is potentially valuable for researchers in gastroenterology and intestinal stem cell biology. However, the manuscript is largely descriptive, and several concerns regarding the experimental design and data interpretation remain:\n1. The cytokine concentrations used appear arbitrary and require justification, ideally supported by preliminary data or references. Furthermore, Table S14 (referenced in line 190), which should contain the concentrations and vendors of the 81 niche factors, is missing. This information is essential for assessing the validity of the findings.\n2. Related to the above, please provide a rationale for the concentrations used in the DSS and irradiation experiments.\n3. The organoid culture media includes Wnt, R-spondin, Noggin, and A83-01, which may interfere with differentiation and block specific cytokine signaling pathways. For example, Noggin inhibits BMP signaling, and A83-01 blocks TGF-β signaling. Please discuss how these factors may influence the experimental outcomes.\n4. Line 309: The hypothesis that IFN-α may have protective effects in inflamed epithelium is intriguing but appears to contradict the known effects of JAK inhibitor treatment in ulcerative colitis. Such a hypothesis should be supported by direct biological experimentation.\n5. Line 151: Is CUCY2C truly unrelated to steroid treatment? Please clarify or correct.\n6. Figure 5B lacks scale bars, which are necessary for proper interpretation of the imaging data."}, "author_response": "Reviewer #1\nThis study presents a valuable single-cell and spatial transcriptomic atlas of human colonic\nepithelial responses to inflammatory stimuli and niche factor perturbations. The manuscript\noffers a comprehensive resource that maps epithelial responses to microenvironmental cues in\nan IBD-relevant organoid model, effectively linking in vitro findings to human disease through\nscRNA-seq and spatial analyses. Non-biased and comprehensive study is entirely valuable data\nsource for proceeding in the related field.\nOne comment is the connection with fetal reversion. or about recent. The authors identify a\nCXCL8⁺AREG⁺ANXA1⁺ population in dissociation-induced organoids (Fig. 1F), which\nresembles previously described revival stem cells (revSCs) as reported in several seminal\nstudies (Yui et al., Cell Stem Cell, 2018; Nusse et al., Cell Stem Cell, 2018; Ayyaz et al., Nature,\n2019; Klementaitė et al., Nature, 2020). However, the potential presence and physiological\nrelevance of such revSC-like populations in IBD patient tissues are not explored. Given that the\nauthors have access to scRNA-seq and CosMx spatial transcriptomic data from ulcerative colitis\npatients, they are in a strong position to assess the expression and localization of canonical\nrevSC markers such as ANXA1, TACSTD2, and CLU in vivo. Incorporating such an analysis\nwould significantly strengthen the interpretation of the regenerative programs described in the\norganoid system.\nIn addition, the gene programs (GPs) extracted in this study are mainly based on in vitro data\ngenerated under a specific matrix condition: all organoids were cultured in 2% Matrigel\nsuspension across injury and perturbation conditions. As shown in Figure 1B, the epithelial cells\nare basically polarized as written in the text. While this uniform setup facilitates comparative\nanalysis, it is important to recognize that matrix composition and mechanical context can\nprofoundly influence epithelial cell states, particularly under regenerative or inflammatory\nconditions. For example, Yui et al. (2018) showed that revSC signatures emerge preferentially\nin collagen-based cultures with the features of depolarization of epithelial cells. Moreover, the\nauthors themselves highlight “cell–matrix adhesion” as a defining component of GP6 (Fig. 4D),\nsuggesting ECM-dependent transcriptional responses.\nI do not think that additional experiments are required. However, I would recommend that the\nauthors consider the following:\n1. Acknowledge the limitations of their matrix system in shaping regenerative cell states. For\nexample, human colon collagen-based organoids have been shown to be enriched in several\nCXCL genes (Kobayashi S. et al., Inflammation and Regeneration, 2022), which may not be\nfully captured in Matrigel-based systems.\nWe acknowledge that our organoid system, like any in vitro model, is a minimalistic system that\nmay not fully recapitulate all aspects of human biology. In particular, it is a good point that the\nchoice of matrix in particular plays a large role in organoid biology. While we show in Extended\nData Fig. 1 that organoids in microwell suspension culture with 2% Matrigel largely resemble\ntraditional culture in Cultrex domes, there are some differences in expression of matrix signaling\nproteins between the two culture systems. In particular, microwell organoids express higher\nlevels of ITGA2, ITGA6, ITGAV, and ITGB1 across 3 experiments (Figure R1). Thus, we have\nincluded text in the Discussion highlighting the limitations of our system in comparison to\ncollagen-based organoids, and have cited the above reference.\nFigure R1: Expression of matrix signaling proteins from scRNA-seq data comparing organoids cultured in 3D Cultrex\ndomes vs. microwell suspension culture with 2% Matrigel.\n2. Discuss the possibility that certain regenerative programs—such as revSC-like states—may\nbe underrepresented or modulated by the current culture context. Alternatively, if such states\nare embedded within specific GPs identified in this study, it would be valuable to relate those\nprograms to the known revSC phenotype and situate the findings within the broader context of\nintestinal inflammation and epithelial regeneration. Even though the perturbation screen is non-\nbiased, the resulting transcriptional programs show substantial overlap with features of\npreviously described revSC states.\nWe found that while some of our gene programs derived from injured organoids contained a few\ngenes that overlap with revSC-like signatures (e.g. ANXA1 and TNFRSF12A in GP19,\nTACSTD2 in GP15, CD44 in GP2, LRIG1 in GP1, LY6D in GP8), none of the GPs we identified\nfully captured the revSC state. Instead, we decided to investigate the revSC signature\nseparately given the presence of the CXCL8⁺AREG⁺ANXA1⁺ population in dissociation-induced\norganoids that seemed to bear relevance to revSCs. Based on gene expression, we found that\ndissociated organoids best elicit the revSC signature, as they upregulate ANXA3, TNFRSF12A,\nCLU, ANXA1, and TACSTD2. However, none of the injury conditions increased expression of all\ngenes associated with revSCs (ANXA1, TACSTD2, CLU, ANXA3, CD44, CCND1, LRIG1, BMI1,\nTERT, PTEN, PPM1D, LY6D, CXADR, TNFRSF12A). Notably, the cytomix injury condition,\nwhich bears the most transcriptional similarity to IBD, does not appear to elicit a revSC-like\nphenotype in human organoids.\nUpon suggestion, we have also explored expression of the revival stem cell signature in both\nhuman patient scRNA-seq and CoxMx spatial transcriptomic data. In the scRNA-seq data, we\nfound that some genes associated with revSC’s are upregulated in inflamed vs. healthy\npatients, although this association seemed higher in the Kong et al. dataset vs. the Smillie et al.\ndataset (Extended Data Figure 4A). In the CosMx dataset, we unfortunately were limited to 8\ngenes that overlapped with the CosMx panel (ANXA1, CCND1, CD44, CLU, LY6D, PTEN,\nTACSTD2, and TNFRSF12A). While we found that this signature as a whole was upregulated in\nUC vs. healthy tissue (Extended Data Figure 4C), the overall expression of the signature was\nlow in the epithelium (Figure R2). Additionally, the signature seemed to be driven by ANXA1,\nCD44, and TNFRSF12A, while other genes showed little change between UC vs. healthy\n(Figure R3).\nFigure R2: Representative CosMx images of UC tissue, highlighting the expression score of the revival stem cell\nsignature.\nFigure R3: Mean expression of individual genes in the revSC signature per sample in UC vs. healthy tissue.\nBased on the potential disease-relevance of this signature, we have included this new data\nfocusing on revival stem cells in Figure 2A and Extended Data Figure 4A and C.\n3. Include key references on revSC biology (Yui et al., 2018; Nusse et al., 2018; Ayyaz et al.,\n2019; Klementaitė et al., 2020 etc., if any others important) to better contextualize the\nregenerative populations described here and enhance the biological relevance of the findings.\nIncorporating these points would improve mainly the Discussion section by providing a more\ncomprehensive and physiologically informed interpretation of the epithelial regenerative\nresponses uncovered in this study.\nThank you for highlighting key references on revSC biology! We had already cited Ayyaz et al.,\n2019 in regard to the CXCL8⁺AREG⁺ANXA1⁺ population in dissociation-induced organoids, but\nhave updated our references to include Yui et al., 2018 and Nusse et al., 2018. We have\nadditionally added more text specifically highlighting the revSC signature when discussing the\ndissociation injury model in the Results, as well as to the Discussion section, particularly in\nregard to the new plots added to Figure 2A and Extended Data Figure 4A and C.\nReviewer #2\nCapeling MM et al. present single-cell transcriptomic data of human colonic organoids subjected\nto various cytokine stimulations. They also integrate their single-cell dataset with spatial\ntranscriptomic data from inflammatory bowel disease (IBD) tissues. This dataset is potentially\nvaluable for researchers in gastroenterology and intestinal stem cell biology. However, the\nmanuscript is largely descriptive, and several concerns regarding the experimental design and\ndata interpretation remain:\n1. The cytokine concentrations used appear arbitrary and require justification, ideally supported\nby preliminary data or references. Furthermore, Table S14 (referenced in line 190), which\nshould contain the concentrations and vendors of the 81 niche factors, is missing. This\ninformation is essential for assessing the validity of the findings.\nThank you for pointing out this oversight, we apologize that we unfortunately missed the upload\nof Table S14 (now Table S13, as we have condensed the previous Table S2 into the Source\nData) upon initial submission. This table includes information on all of the niche factors included\nin the screen, including the vendor, catalog number, and concentration used in our study.\nThe cytokine concentrations were selected to be within or above the EC listed from the vendor\n50\nproduct information in order to increase the chances of observing a transcriptional effect, which\nfor most compounds was 100ng/mL. For our pro-inflammatory cytokines TNFɑ and IFNɣ, we\nhave observed evidence of toxicity and cell death in preliminary experiments at 100ng/mL, so\nwe reduced the concentrations of these molecules down to 50ng/mL. Additionally, we increased\nthe concentrations for some factors (ANGPTL1, ANGPTL4, AREG, BDNF, BMP8B, CT-1,\nCTSG, F2, FGF1, FGF18, FGF7, FGF9, NT-4, POSTN, PPY, TWEAK, VIP) up to 250ng/mL\nbased on preliminary data in our lab showing that these compounds did not elicit strong\ntranscriptional effects at 100ng/mL. While we understand that EC values are cell type\n50\ndependent, it unfortunately is not tractable to do dose titrations with a scRNA-seq readout, nor is\nthere a standard assay to measure pathway activation in an efficient manner for all ligands\nused. We have added language describing this justification to the Methods section to increase\nclarity.\n2. Related to the above, please provide a rationale for the concentrations used in the DSS and\nirradiation experiments.\nFor the DSS experiments, we selected the concentrations of 0.1% and 0.5% based on\npreliminary experiments in which we found that higher concentrations of DSS (1% or greater)\nlead to a significant decrease in the expression of the stem cell marker LGR5 after 24 hours as\nassessed by qRT-PCR, and that LGR5 expression was not able to be recovered after 3 days.\nAdditionally, 1-2% DSS caused significant organoid death. However, with 0.1% DSS, LGR5\nexpression was significantly reduced after 24 hours but recovered 3 days after treatment, and\norganoids remained largely viable. Thus, we decided to move forward with testing our ‘low dose’\nDSS as 0.1% such that organoids would exhibit a recoverable defect in stem cell number, and a\n‘high dose’ as 0.5% to remain below the toxic dose of 1%.\nFor the irradiation experiments, we selected concentrations of 6Gy and 12Gy based off of a\nstudy assessing the impact of ɣ-irradiation in mouse intestinal organoids (Montenegro-Miranda\net al, 2020) and based on our preliminary experiments. In our pilot data, we found that 6Gy and\n10Gy irradiation caused a decrease in expression of the proliferation marker MKI67 compared\nto controls or 2Gy irradiation as assessed by qRT-PCR, suggesting that doses of 6Gy or higher\nelicit an injury response in organoids. However, since we did not see a strong transcriptional\nresponse to 6Gy irradiation, we decided to further increase the concentration of our ‘high’ dose\nirradiation to 12Gy.\nFor the cytomix experiment, we performed scRNA-seq on organoids treated with varying\nconcentrations of TNFɑ, IFNɣ, and IL1β and assessed the expression of gene programs IBD1\nand IBD2 in response to treatment. Here, we found that expression of these IBD-related GPs is\ndose dependent, and selected a low and high dose based on differential expression of these\nGPs.\nTo add clarity surrounding the injury doses selected for inclusion, we have added qRT-PCR and\nscRNA-seq data from our pilot experiments to Figure S2A and have added text to the Methods\nclarifying that doses were selected on preliminary dose-response studies.\n3. The organoid culture media includes Wnt, R-spondin, Noggin, and A83-01, which may\ninterfere with differentiation and block specific cytokine signaling pathways. For example,\nNoggin inhibits BMP signaling, and A83-01 blocks TGF-β signaling. Please discuss how these\nfactors may influence the experimental outcomes.\nOur organoid culture media includes factors that promote organoid growth and maintenance,\nwhich are commonly used in the field. For a screen involving many growth factors, we decided\nto maintain our organoid base media with all such factors rather than withdrawing them to\nmaintain organoid health, especially since our screen already involved the addition of pro-\ninflammatory cytokines that cause stress and cell stress and cell death. Nonetheless, it is a valid point that the\ninclusion of Noggin and A83-01 may interfere with the addition of recombinant BMP and TGF-β\nfamily ligands. However, as our secreted factors were added to the culture media at the upper\nend of active/published concentrations, we believe that the doses we used are high enough to\ncounteract inhibitors already in the media. Noggin has been shown to primarily inhibit BMP2 and\nBMP4 with weaker effects on BMP6, but BMP2 and BMP4 were some of the perturbations with\nthe strongest transcriptional effect and highest number of DEGs, comparable to BMP6 which\nshould have minimal inhibition from Noggin, and greater than other perturbations like IL22 or\nHBEGF which do not face inhibition from any media components. Similarly, we see a high\ntranscriptional response and number of DEGs for TGF-β1,2, and 3, suggesting that we are\noperating at a dose high enough to counteract inhibition from A83-01.\nBased on these results, we believe that the high concentrations of secreted factors used should\nbe enough to counteract inhibitors in the media. For transparency, we have added a statement\ncommenting on the presence of such inhibitors in the media to the Discussion section. We also\ndiscuss that if any of these above perturbations are of significant interest to the reader, follow up\nexperiments are warranted to confirm DEGs are conserved in presence or absence of media\nmorphogens.\n4. Line 309: The hypothesis that IFN-α may have protective effects in inflamed epithelium is\nintriguing but appears to contradict the known effects of JAK inhibitor treatment in ulcerative\ncolitis. Such a hypothesis should be supported by direct biological experimentation.\nWe agree that the overall transcriptional impact of IFN-α appears to be more pro-inflammatory\nthan protective. In Figure 5D, we show that IFN-α elicits upregulation of GP IBD1 and\ndownregulation of IBD2, matching the expression pattern of cytomix treatment and IBD. In order\nto further investigate the impact of IFN-α vs. JAK inhibition on colonic organoids, we conducted\na follow-up experiment where we treated colon organoids from 3 donors, either control or\ncytomix-treated, with IFN-α and the JAK inhibitor Tofacitinib. In this experiment, we again found\nthat IFN-α upregulates GP IBD1 and downregulates GP IBD2, to an even stronger extent than\ncytomix treatment alone, whereas JAK inhibition clearly shows downregulation of GP IBD1 and\nupregulation of IBD2 (Figure R4). While previous studies have noted that IFN-α may exhibit\ncontext-dependent roles to promote inflammation or homeostasis (Rauch, I. et al. Eur. J.\nImmunol. 2014), we did not find enough evidence to fully support a protective role for IFN-α,\nespecially since the expression pattern of IBD-related gene programs IBD1 and IBD2 align with\nthe direction of inflammation in two experiments.\nFigure R4: Expression of gene programs IBD1 and IBD2 in organoids treated with cytomix, IFN-α and Tofacitinib in a\nfollow up experiment\nBased on these results, we have decided to remove language suggesting a protective role for\nIFN-α since such a conclusion may be misleading based on the overall transcriptional profile of\nIFN-α. Instead, we have decided to highlight another organoid-derived gene program in Figure\n5, GP12, which shows stronger effect sizes for multiple perturbations and allows us to highlight\nways in which we can infer the directionality of transcriptional responses based on our control\nperturbations.\nWe found that GP12 (related to epithelial structure maintenance and negative regulation of\nanoikis) showed relevance to IBD, as this signature was upregulated in IBD vs. healthy\nepithelium in both scRNA-seq and spatial transcriptomic data. However, this signature was\nactually downregulated with cytomix/pro-inflammatory cytokines in our organoid data. Since we\nknow that cytokines such as TNFɑ and IFNɣ drive inflammation, this suggests that upregulation\nof GP12 may be protective in the context of human tissue. We have adjusted the text and\nFigure 5 accordingly to eliminate confusion, as further investigation to piece apart the potential\ninflammatory vs. protective roles of IFN-α is outside the scope of the current work.\n5. Line 151: Is CUCY2C truly unrelated to steroid treatment? Please clarify or correct.\nLine in question :\nBased on our gene ontology analysis using the \"GO_Biological_Process_2023\" database\nimplemented in the enrichR R package, GUCY2C is indeed annotated as being associated with\n“Steroid Metabolic Process (GO:0008202)”. We have added the GO terms associated with\ngenes upregulated by BMP2/6 in cytomix and EGF/HBEGF in DSS to Table S5 for reference.\nHowever, upon closer inspection, we agree that GUCY2C is not the most representative gene\nfor that GO term. Indeed, BMP2 and 6 also upregulate more canonical genes associated with\nsteroid metabolism, including HSD11B2, AKR1C2, and AKR1C3. For clarity and accuracy, we\nhave updated the genes described in the text.\n6. Figure 5B lacks scale bars, which are necessary for proper interpretation of the imaging data.\nThank you for pointing out this oversight! We have added the proper scale bars to our CosMx\nimages."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have provided a comprehensive and thoughtful rebuttal, and the revised manuscript has improved substantially.\nI appreciate the extensive additional analyses they performed—particularly the examination of revSC-associated genes in both scRNA-seq datasets (Kong and Smillie) and CosMx spatial transcriptomics, as well as the expanded discussion of matrix-dependent effects on epithelial states. These additions directly and thoroughly address the main points raised in my initial review.\nRegarding the revival stem cell (revSC)–like biology, the authors convincingly demonstrate that while certain injury conditions or dissociation-induced states upregulate subsets of canonical markers, no single condition fully captures the classical revSC signature. Their in vivo analyses further indicate that the expression of these markers is detectable but overall limited within UC patient epithelium, driven mainly by a subset of genes (e.g., ANXA1, CD44, TNFRSF12A). This modest overlap may reasonably reflect biological heterogeneity across patient cohorts, inflammatory states, and tissue sampling contexts, and I appreciate that the authors now discuss these considerations explicitly.\nI also commend the authors for incorporating key references in revSC biology and for expanding the Discussion to more clearly contextualize their findings within the broader literature on regenerative epithelial states. Their acknowledgement of matrix-related limitations—along with the additional data comparing 3D domes and 2% Matrigel microwell cultures—provides appropriate transparency and strengthens the interpretation of gene programs derived from the perturbation screen.\nOverall, the authors have addressed all my comments with rigor and clarity.", "reviewer_2": "The authors have addressed concerns I raised before."}, "author_response": ""}]

Decision letter