Organization of the human intestine at single-cell resolution

Abstract

The intestine is a complex organ that promotes digestion, extracts nutrients, participates in immune surveillance, maintains critical symbiotic relationships with microbiota and affects overall health¹. The intesting has a length of over nine metres, along which there are differences in structure and function². The localization of individual cell types, cell type development trajectories and detailed cell transcriptional programs probably drive these differences in function. Here, to better understand these differences, we evaluated the organization of single cells using multiplexed imaging and single-nucleus RNA and open chromatin assays across eight different intestinal sites from nine donors. Through systematic analyses, we find cell compositions that differ substantially across regions of the intestine and demonstrate the complexity of epithelial subtypes, and find that the same cell types are organized into distinct neighbourhoods and communities, highlighting distinct immunological niches that are present in the intestine. We also map gene regulatory differences in these cells that are suggestive of a regulatory differentiation cascade, and associate intestinal disease heritability with specific cell types. These results describe the complexity of the cell composition, regulation and organization for this organ, and serve as an important reference map for understanding human biology and disease.

Fig. 1. CODEX multiplexed imaging of eight regions from the small intestine and colon to create a single-cell map of the healthy human intestine.

a–c, Cell type percentages from CODEX data averaged across eight donors. Cell types are normalized to the stromal (a), immune (b) and epithelial (c) compartments. Statistical analysis was performed using two-sided t-tests comparing the difference in cell type percentage between the small bowel (SB) and the colon (CL); *P < 0.05. ICC, interstitial cells of Cajal; NK, natural killer cells; TA, transit-amplifying cells. d, The percentage of M1 macrophages within the small bowel and colon for all donors plotted versus donor BMI (Pearson correlation r = 0.86). e, Cell type maps of the mid-jejunum from representative individuals (n = 8 donors) with high or low BMI with M1 macrophages (black) highlighted among stromal (light grey) and epithelial (grey) cell types also shown. Scale bar, 250 µm. f, Cell type percentages for endothelial and CD8⁺ T cells compared for donors with or without a history of hypertension. Statistical analysis was performed using two-sided t-tests; *P = 0.038, ***P = 0.00013. n = 3–5 donors. g, Quantification of the same-cell density measured as an average distance of its five nearest same-cell neighbours divided by the maximal possible same-cell distance within the tissue. n = 64 tissue sections. h, The percentage of macrophage subsets across major intestinal compartments.

Fig. 2. Multicellular neighbourhood analysis of the intestine.

a, Twenty unique intestinal multicellular neighbourhoods were defined by enriched cell types as compared to the overall percentage of cell types in the samples. b, An example in which neighbourhoods mapped back to the tissue show overall tissue structures. Scale bar, 0.5 mm. c, The percentage of Neuroendocrine-Enriched neighbourhood of all of the neighbourhoods as determined by individually characterizing cellular neighbourhoods by region. n = 8 donors. The box plots show the median (centre line), 25th to 75th percentile (box limits), minimum and maximum values (whiskers), and outliers (points outside 1.5× the interquartile range). d,e, Quantification of the same-cell density for neuroendocrine cells compared across the small bowel and colon (n = 32 tissue sections) (d) or the epithelial neighbourhoods as determined by individually characterizing cellular neighbourhoods by region (n = 64 tissue sections) (e). Avg., average; max., maximum. f, A subset of epithelial neighbourhoods mapped back to a representative magnified region (n = 8 donors) of the mucosa of a transverse colon section. Scale bar, 250 µm. g,h, The approach to calculate concentric increasing neighbourhoods around a Paneth cell (g) to generate cellular neighbourhoods for Paneth cells at increasing radii (h). i,j, Schematic (i) and CODEX fluorescence data illustrating a representative (1 of 32 sections from 8 donors) magnified portion of the proximal jejunum depicting colocalization of Paneth cells (DEFA5, green) and CD8⁺ T cells (CD8, cyan) and CD4⁺ T cells (CD4, yellow) in the intestinal crypt environment (j). Scale bar, 50 µm.

Fig. 3. Multilevel hierarchical structural description of the small intestine and colon.

a–d, Representation of multiple levels of hierarchical description: cell type (a), multicellular neighbourhood (b), community (based on clustering windows of cell neighbourhoods) (c) and tissue units (based on clustering communities) (d) comparing the small bowel with the colon for two representative tissue sections (from a total of 64 sections from 8 donors). Scale bar, 1 mm. e, Tissue hierarchy graph of the multilevel network of the tissue comprised of the different structures. Shapes correspond to structural level (cell type, neighbourhoods, communities, tissue units); colours represent individual categories as indicated in a–d; the size of shapes represents the percentage of tissue; and the size of connected lines represents the overall contribution to the next level of the structure when moving down the graph in increasing tissue structural hierarchy. The black rectangles highlight a single trajectory highlighted within this Article. The red bracket indicates separation of stromal tissue units from the mucosal tissue units. f, Magnified mucosal area of a colon community map shown within c. Scale bar, 100 µm. g, The spatial-context maps of the colon highlighting relationships of communities across the entire sample. This structure is defined by the number of unique communities required to make up at least 85% in a given window. The circles represent the number of cells represented by a given structure. The green rectangle highlights a structure discussed in this Article and maps this structure back to g. The colours are as indicated in c. h, The cell type percentage of immune cells shown for each community ordered in relative order of general increasing proximity to the lumen on the basis of community spatial-context analysis.

Fig. 4. Single-cell atlas of gene expression and chromatin accessibility in the human intestine.

a, Sections of the intestine analysed by separate snRNA-seq and snATAC-seq experiments or multiome experiments. b,c, UMAP representation of all snRNA stromal (b) and immune (c) cells coloured by cell type. DC, dendritic cells; fib., fibroblasts; ILC, innate lymphoid cells; myofib., myofibroblasts; SM, smooth muscle cells. d, UMAP representation of snRNA epithelial cells in the four primary regions of the intestine. Jejunum includes both proximal- and mid-jejunum samples. Colon includes samples from the ascending, transverse, descending and sigmoid colon. e, Expression of INSL5 and the INSL5 receptor, RXFP4, in different cell types in different regions of the intestine. RXFP4 was expressed in less than 2.5% of all epithelial cell types that were not included in the dot plot. C, colon; D, duodenum; I, ileum; J, jejunum; enterochrom., enterochromaffin cells; enteroendo., enteroendocrine cells. f, Beeswarm plot showing the log-transformed fold change between the small intestine and colon for groups of nearest-neighbour cells from different cell type clusters. Significant changes are indicated in red and blue. Lymph. endo., lymphatic endothelial cells. g, Subclustering of specialized secretory cells in d coloured by cell type. h, The expression of secretory genes in specialized secretory cells defined in g. i, The percentage of MUC6⁺ enterocytes among all cell types for the four donors imaged using MUC6 antibodies. j, Representative CODEX fluorescence image of the duodenum (6 out of 57 markers overlaid) (left). Hoechst (nuclei), MUC6, MUC1 (also found in gland areas), cytokeratin (pan-epithelial), α-SMA (muscle) and vimentin (stromal) staining is shown. Right, the magnified area highlights the gland just below the mucosa in the submucosa. This experiment was independently repeated four times. Scale bars, 500 µm (left) and 50 µm (right). Source Data

Fig. 5. Integration of CODEX multiplexed imaging, snRNA-seq and snATAC–seq reveals genes associated with cellular neighbourhoods and regulatory TFs in the human intestine.

a, snRNA-seq and CODEX integration using MaxFuse for both the small bowel and colon. Integrated cells are coloured by modality of origin, cell types or RNA/protein expression levels. MUC2 expression (protein and RNA) is shown on the integrated plots, as well as the top 10 DEGs of the Inner Follicle neighbourhood as determined by CODEX. Min., minimum. CN, cellular neighbourhood. b, TFs of which the integrated gene expression correlated with their motif activity in one region of the intestine. Row z-scores of ChromVar deviation scores are shown on the left and row z-scores of integrated TF expression are shown on the right. Ent., enterocytes. abs., absorptive. Source Data

Fig. 6. Regulation of differentiation in the human intestine.

a, UMAP projections depicting the cells in the four primary regions of the intestine (duodenum, jejunum, ileum and colon), labelled by cell type (left) and differentiation pseudotime (right). b,c, Variable peaks (b) and genes (c) identified along the absorptive differentiation trajectories. The rows represent the z-scores of accessibility for each peak or expression for each gene. The columns represent the position in pseudotime from the start to the end for each section of the intestine. Peaks and genes were k-means-clustered and the clusters were labelled on the basis of the dominant time and location where they are most accessible/expressed. d, Integrated gene expression of TFs of which the expression is correlated with ChromVar motif activity along the differentiation trajectory. e, Accessibility at peaks correlated with the expression of ETV6 along the differentiation trajectory in each region (left). Each peak is normalized to the maximum accessibility along any of the trajectories. Right, integrated gene expression of ETV6 along the differentiation trajectory in each region is plotted on the right. Norm., normalized. f, Linkage-disequilibrium score regression to identify the enrichment of GWAS SNPs in cell-type-specific marker peaks. Unadjusted coefficient P values computed from linkage-disequilibrium score regression are plotted in the heat map. Significance is indicated by an asterisk, as determined by a Bonferroni-corrected coefficient P value of <0.05. P values for determining significance were adjusted for the number of cell classes tested. Source Data

Extended Data Fig. 1. CODEX mulitplexed imaging across the healthy human intestine reveals changes in cell composition and organization.

A) Schematic for how CODEX multiplexed imaging was performed on arrays of 4 different sections of either colon and small intestine from the same donor simultaneously. Image processing steps done to extract single-cell spatial data. B) An example CODEX fluorescent image of one region of the small bowel (SB) (1 of 64 tissue sections) for CODEX with 6/54 markers shown for one donor (scale bar = 1 mm and magnified insert = 100 µm) with C) accompanying cell type map following cell segmentation and unsupervised clustering. D-E) Stromal cell type percentages either as a percent of D) All stromal cells, or E) all cells restricted to the Muscularis Externa tissue unit. F) Immune cell type or G) epithelial cell type percentages either as a percent of all cells restricted to the Mucosa tissue unit. H) Percentage of CD57+ Enterocyte cells of all cell types across different areas samples from small intestine to colon. (for D-H: * p value< 0.05, ** p value< 0.01, *** p value < 0.001 by two-sided T test, n=8 donors). (All boxplots in figures are plotted as minimum, 25 percentile, median, 75 percentile, maximum, and outliers as points outside 1.5 the interquartile range). I) Cell map of a representative section (one of 8 donors) of the Duodenum that shows CD57+ Enterocyte presence in glands in the Submucosa where Enterocytes and TA cells are shown in dark grey and Smooth muscle cells in light grey (other cell types not shown) (scale bar = 500 µm). J-K) Quantification of the same-cell density that is measured as an average distance of its 5 nearest same-cell neighbours normalized by the maximal possible same-cell distance within the tissue (n = 64 tissue sections) for J) stromal and K) epithelial cell types. L) A representative cell type map (one of 64 tissue sections from 8 donors) with only plasma cells, CD8+ T cells, and M2 Macrophages shown (scale bar = 500 µm).

Extended Data Fig. 2. Multicellular neighbourhood analysis of CODEX imaging data indicates conserved cellular structures across the intestine.

A) Neighbourhood analysis was done by taking a window across cell type maps and vectorizing the number of cell types in each window, clustering windows, and assigning clusters as cellular neighbourhoods of the intestine. B-D) Neighbourhood percentages from CODEX data averaged normalized by B) stromal, C) immune, and D) epithelial compartments. Asterix indicates p-value less than 0.05 difference in cell type percentage from the small bowel (SB) to the colon (CL) by two-sided T test. E) Stromal multicellular neighbourhood percentages either as a percent of all neighbourhoods restricted to the Muscularis Externa tissue unit (* p value< 0.05, *** p value < 0.001, n = 8 donors, by two-sided T test,). F) Quantification of the same-cell density for just smooth muscle cells within different smooth muscle multicellular neighbourhoods (x axis) (n = 32 tissue sections). G-H) Cell type maps for a region of the small intestine (one of 64 tissue sections imaged from 8 donors) with G) all cell types plotted for the whole tissue (scale bar = 500 µm), H) cells contained within the plasma cell neighbourhood (scale bar = 500 µm), and a magnified area of denoted by rectangle showing subset of cell types (scale bar = 50 µm). I) CODEX fluorescent imaging with subset of fluorescent markers overlaid for the same tissue as G (Hoechst=Blue, CD4=Green, CD68=magenta, CD38=yellow, CD206=cyan, CD138=grey), (scale bar = 500 µm with magnified insert scale bar = 100 µm).

Extended Data Fig. 3. Change in multicellular neighbourhoods across the intestine.

A) Immune neighbourhood percentages as a percent of immune neighbourhoods. B) Neighbourhood percentages for Microvasculature and CD8+ T cell IEL neighbourhoods compared for donors with or without a history of hypertension (Microvasculature p value = 0.0065<, CD8+ T Enriched IEL p value = 0.0017 by two-sided T test, n = 3-5 donors). C) Difference in composition in neighbourhood by cell type for all neighbourhoods based on subtracting the log2 fold enrichment of each cell type found within that neighbourhood compared to average percentages in the tissue in SB from CL. Neighbourhoods and cell types are ordered by summing the absolute value of all rows and columns to denote conservation of a neighbourhood. D) 22 unique intestinal multicellular neighbourhoods (y axis of heatmap) were defined by enriched cell types (x axis of heatmap) as compared to overall percentage of cell types in the samples with 2 unique neighbourhoods not identified with overall neighbourhood analysis. E) Epithelial neighbourhood percentages as a percent of epithelial neighbourhoods from multicellular analysis performed on each individual region of the intestine separately. F) Immune neighbourhood percentages as a percent of immune neighbourhoods from multicellular analysis performed on each individual region of the intestine separately. (* p value< 0.05, ** p value< 0.01, *** p value < 0.001 by two-sided T test, n = 8 donors).

Extended Data Fig. 4. Multi-neighbourhood community analysis of the intestine.

A) Community analysis was done by taking a window across neighbourhood maps and vectorizing the number of each neighbourhood type in each window, clustering windows, and assigning clusters as multi-neighbourhood communities of the intestine. 10 unique intestinal multi-neighbourhood communities (y axis of heatmap) were defined by enriched neighbourhood types (x axis of heatmap) as compared to overall percentage of neighbourhood types in the samples. B) Quantification of neighbourhood types across each section of the intestine (colour legend within panel A). C) Community percentages as a percent of all communities for small intestine and colon (* p value< 0.05, ** p value< 0.01, *** p value < 0.001, n = 8 donors). D) Concentric multi-neighbourhood analysis surrounding only neighbourhoods labelled as Paneth Cell Enriched neighbourhoods, with number of nearest neighbours for a given row in the heatmap.

Extended Data Fig. 5. Multicellular neighbourhood and community interactions across the intestine.

A) Schematic of a representative community map and corresponding tissue community network graph (one of 64 tissues). B-C) Community-community motifs that are significantly enriched in both the small intestine and colon as compared to a null distribution of motif instances created from random permutation of tissue graph labels, where B) shows only those motifs that interact with the Adaptive Immune Enriched community and C) shows all other shared motifs between the SB and CL. Motifs indicated by shape and colour indicate those motifs that have significant p value versus those that are indicated with an x in the graph; p values were Bonferroni corrected by multiplying by twice the number of tests conducted in each comparison group. Colour legend is also the same as panel B. D) Representative neighbourhood map (one of 64 tissue sections from 8 donors) with (scale bar = 500 µm) E) Region magnified as in the main figure of the mucosal area of a colon community map, but this time with the multicellular neighbourhoods coloured (see panel D for legend) (scale bar = 100 µm). F) Spatial context maps of the CL highlighting relationships of multicellular neighbourhoods across just the neighbourhoods found within the tissue unit Mucosa. This structure is defined by the number of unique neighbourhoods required to make up at least 85% in a given window. Circles represent the number of cells represented by a given structure. Red rectangle highlights a structure discussed in the manuscript and maps this structure back to panel K. Colour legend is also the same as panel D. G-H) Cell type percentage of G) epithelial and H) stromal cells shown for each community ordered in relative order of general increasing proximity to the lumen based on community spatial context analysis.

Extended Data Fig. 6. Quality control and clustering of single-nucleus data.

(a,b) Violin plots of TSS enrichment (a) and RNA counts/cell (b) for different samples included in the study. Samples are coloured by the location from which they were obtained. (c) UMAP projection of all snRNA cells coloured by location. (d) Dotplot representation of expression of immune marker genes by immune cell types. (e) UMAP projection of scRNA immune cells coloured by donor. (f) Dotplot representation of expression of stromal marker genes by stromal cell types. g) Sub-clustering of enteroendocrine cells from all regions of the intestine. h) Dotplot representation of the expression of subtype specific enteroendocrine and enterochromaffin marker genes in different enteroendocrine cell types in our datasets. I) Sub-clustering of specialized secretory cells coloured by sample.

Extended Data Fig. 7. Differential cell type abundance.

a–c) Beeswarm plot showing the log-fold change between the three main regions of the small intestine and colon for groups of nearest neighbour cells from different cell type clusters in the stromal (a), immune (b), and epithelial (c) compartments computed with Milo. Significant changes are indicated in colour. d) Boxplots comparing the fraction of all cells in each sample composed of each cell type for samples from the colon, ileum, jejunum, and duodenum. e) Log2FC in abundance of each cell type between the regions listed on the y-axis as estimated with scCODA. Only significant results at an FDR of 0.05 are shown, with all nonsignificant differences plotted as white.

Extended Data Fig. 8. Pairing of CODEX multiplexed imaging and snRNA-seq for cell-cell colocalization analysis.

a) Expression of B3GAT1 and four CEACAM transcripts plotted on the UMAP manifold of epithelial cells from the duodenum, jejunum, ileum, and colon. b) Dotplot representation of the expression of B3GAT1 and four CEACAM transcripts by different epithelial cell types in different regions of the intestine. c) Large colon (CL) and small bowel (SB) show differences in cell-cell co-localization patterns; annotated cell-pairs are more colocalized in the colon compared to the small bowel (Student’s T test, two-sided, corrected for multiple hypothesis testing with the Benjamini Hochberg procedure). d) SEMA4D ligand expression in plasma cells and MET receptor gene expression in TA2 cells, showing higher expression in colon than small bowel (one-sided Wilcoxon Rank-Sum Test). e) Differences in pairwise cell-type colocalization patterns across tissue locations (n = 3 samples for each location).

Extended Data Fig. 9. Integration of snRNA-seq and CODEX datasets.

a) Average MET expression for all TA cells from a given donor and average SEMA4D expression for all plasma cells from a given donor. b) PLXNA2, RASA1, and MAP2K5 expression in TA2 cells in large colon (CL) and small bowel (SB) (one-sided Wilcoxon Rank-Sum Test). c) Representative image of a donor’s transverse colon with plasma cells (red), TA cells (blue), and other cell types (grey) highlighted with also a magnified area indicated with rectangle. (left scale bar = 500 µm; right scale bar = 100 µm) d) Example of receptor ligands expressed at higher levels in the colon than the small intestine (SI) for the FN1 by myofibroblast and PLAUR by enterocyte. e) Colocalization quotient (CLQ) for all cell types found within the follicle community. f) Heatmap of differentially expressed genes (from MaxFuse matched snRNA-seq cells) among individual CODEX cells, grouped based on previously determined cellular neighbourhoods from CODEX analysis. g) Gene Ontology Enrichment analysis for cellular neighbourhoods “inner folliclle” and “outer folliclle”, based on the differentially expressed genes shown in f.

Extended Data Fig. 10. Regulatory TFs in the intestine.

a) Hypergeometric p-values of selected TF motifs in differential peaks between colon stem cells and stem cells from other regions of the intestine. Colour represents the log10 adjusted p-value computed with ArchR. b) Violin plots of motif deviation scores for FOXA3 for goblet cells, immature goblet cells, stem cells, and enterocytes in different regions of the intestine. c) Dotplot representation of expression of different transcription factors in goblet cells, immature goblet cells, stem cells, and enterocytes in different regions of the intestine. d) TF motif footprints for POU2F3 in proximal jejunum tuft cells and enterocytes. e) TF motif footprints for ATOH1 in duodenum goblet cells and enterocytes. Error bands in d and e represent the standard deviation. f) Overlap of epithelial regulators identified with the multiome data and separate snRNA and snATAC datasets.(g, h) Heatmap representation of transcription factors whose integrated gene expression was correlated with their motif activity in one region of the intestine for immune (g) and stromal (h) cell types. Row z-scores of ChromVar deviation scores are shown on the left and row z-scores of integrated TF expression are shown on the right. i, j) Overlap of immune (i) and stromal (j) regulators identified with the multiome data and separate snRNA and snATAC datasets. k) Hypergeometric p-values of TF motifs enriched in the clusters of peaks identified in b. l) Enrichment of KEGG pathways in the clusters of genes identified in c. Uncorrected p-values as determined by kegga are plotted. (m–o) Integrated gene expression of MTTP (m) and SCNN1B (n) and LGR5 (o) along the differentiation trajectory. (p) Accessibility at peaks correlated with the expression of TMPRSS15 along the differentiation trajectory in each region is plotted on the left. Each peak is normalized to the maximum accessibility along any of the trajectories. Integrated gene expression of TMPRSS15 along the differentiation trajectory in each region is plotted on the right.