Single-Cell Analysis Reveals Unexpected Cellular Changes and Transposon Expression Signatures in the Colonic Epithelium of Treatment-Naïve Adult Crohn's Disease Patients

Abstract

Background & aims: The intestinal barrier comprises a monolayer of specialized intestinal epithelial cells (IECs) that are critical in maintaining mucosal homeostasis. Dysfunction within various IEC fractions can alter intestinal permeability in a genetically susceptible host, resulting in a chronic and debilitating condition known as Crohn's disease (CD). Defining the molecular changes in each IEC type in CD will contribute to an improved understanding of the pathogenic processes and the identification of cell type-specific therapeutic targets. We performed, at single-cell resolution, a direct comparison of the colonic epithelial cellular and molecular landscape between treatment-naïve adult CD and non-inflammatory bowel disease control patients.

Methods: Colonic epithelial-enriched, single-cell sequencing from treatment-naïve adult CD and non-inflammatory bowel disease patients was investigated to identify disease-induced differences in IEC types.

Results: Our analysis showed that in CD patients there is a significant skew in the colonic epithelial cellular distribution away from canonical LGR5+ stem cells, located at the crypt bottom, and toward one specific subtype of mature colonocytes, located at the crypt top. Further analysis showed unique changes to gene expression programs in every major cell type, including a previously undescribed suppression in CD of most enteroendocrine driver genes as well as L-cell markers including GCG. We also dissect an incompletely understood SPIB+ cell cluster, revealing at least 4 subclusters that likely represent different stages of a maturational trajectory. One of these SPIB+ subclusters expresses crypt-top colonocyte markers and is up-regulated significantly in CD, whereas another subcluster strongly expresses and stains positive for lysozyme (albeit no other canonical Paneth cell marker), which surprisingly is greatly reduced in expression in CD. In addition, we also discovered transposable element markers of colonic epithelial cell types as well as transposable element families that are altered significantly in CD in a cell type-specific manner. Finally, through integration with data from genome-wide association studies, we show that genes implicated in CD risk show heretofore unknown cell type-specific patterns of aberrant expression in CD, providing unprecedented insight into the potential biological functions of these genes.

Conclusions: Single-cell analysis shows a number of unexpected cellular and molecular features, including transposable element expression signatures, in the colonic epithelium of treatment-naïve adult CD.

Keywords: BEST4; Colonocyte; Crohn’s Disease; Epithelium; Gene Expression; Genome-Wide Association Study; ISC; LGR5; SPIB; Single-Cell; Transposable Element.

Graphical abstract

Figure 1

Clustering and identification of contaminating immune cells. (A) Uniform Manifold Approximation and Projection (UMAP) of initial clusters (left) and clusters following cell type annotation using known cell type markers and highly enriched genes (right). (B) Violin plots showing normalized expression of known markers of different immune cell types on the y-axis and the immune clusters on the x-axis. IEL, intraepithelial lymphocytes; NK, natural killer cells; Pan-B, All B cells; Pan-T, All T cells; Th, T helper cells; Treg, Regulatory T cells.

Figure 2

Single-cell landscape in treatment-naïve Crohn’s disease and NIBD patients. (A) Uniform Manifold Approximation and Projection (UMAP) of epithelial cell clusters following assignment of cell type. (B) Normalized expression heatmap of genes known to be markers of various epithelial cell populations or found to be highly enriched in a cluster (rows) across epithelial clusters (columns) confirms cell type assignment. Normalized gene expression is scaled by row. (C) UMAP overlain with the normalized expression of cell type markers further confirms cluster assignment. MUC2, goblet; LGR5, stem; CA1, CA1+ colonocytes; CEACAM7, CEACAM7+ colonocytes; SPIB, SPIB+ cells; CENPA, G2–M–G1 TA; PCNA, S-phase TA; CHGA, EEC. (D) Normalized expression heatmap of 3 clusters with a colonocytic signature (CA1+ late colonocytes, CEACAM7+ colonocytes, and SPIB+ cells) show expression of cluster-specific markers. Normalized expression is scaled by row. (E) Crypt–axis scores (low near crypt bottom, high near crypt top) of cells shows the expected location of clusters along the crypt axis. Clusters are arranged on the x-axis by mean crypt–axis score (see the Methods section for more detail). The size of the dot corresponds to the level of expression of the known marker of the crypt top: SELENOP. A diagram of colonic crypt is shown (right). (F) UMAP overlain with the RNA velocity. RNA velocity is the ratio of spliced transcripts/unspliced transcripts. Scale is normalized to 1. TA, transit-amplifying.

Figure 3

CD causes shifts in the colonic epithelial landscape. (A) Crypt–axis score density for NIBD and CD cells separately shows a shift toward the colonic crypt top in CD. The significance of the shift was determined using the Kolmogorov–Smirnov test (∗∗∗P < .001). (B) Mean cell abundances across NIBD and CD samples show a significant increase in CA1+ colonocytes in CD. Dots show abundance values of individual samples. (C) Mean cell abundances for either CA1+ early or late colonocyte clusters across NIBD and CD samples shows a significant increase only in the CA1+ late colonocytes. Dots show abundances of individual samples. Abundances of other clusters are not shown. (D) Mean ratio of absorptive to secretory cells across NIBD and CD samples. (E) Mean ratio of EEC to mature goblet cells across NIBD and CD samples shows a significant shift in the secretory lineage toward mature goblet cells in CD. (F) Mean ratio of immature goblet to mature goblet across NIBD and CD samples show a shift toward mature goblet cells in CD. If not indicated, P values are calculated by a Student t test. ∗P < .05; ∗∗P < .01.

Figure 4

Effect of CD within cell types. (A) Number of significant differentially expressed genes in CD vs NIBD. Clusters are shown on the x-axis followed by the number of cells within the cluster in parenthetical notation and the number of differentially expressed genes is shown on the y-axis. (B) Differentially expressed genes for 3 clusters: CA1+ late colonocyte, CEACAM7+ colonocyte, and mature goblet. Log₂ fold change is shown on the x-axis and the adjusted P value is shown on the y-axis. Mean expression of (C) SLC26A2 and CA2, (D) AQP8, (E) CD177 and LYPD8, and (F) CLDN4 across clusters for NIBD and CD samples shows the specificity of expression and differential expression within clusters. P values were calculated by the Wilcoxon rank-sum test. ∗∗∗Adjusted P < .001.

Figure 5

CD disrupts EEC homeostasis. (A) Mean normalized expression across 10 transcription factors involved in EEC maturation in NIBD and CD samples. (B) Mean normalized expression across 7 EEC hormones in NIBD and CD samples. P values were calculated using the Wilcoxon rank-sum test. ∗P < .05, ∗∗P < .01, and ∗∗∗P < .001.

Figure 6

Transposable elements define cell types and are dysregulated in CD. (A) Uniform Manifold Approximation and Projection (UMAP) of epithelial cell clustering after TE inclusion maintains structure. Mean expression of (B) MER11A and (D) MER11C across clusters for NIBD and CD samples shows the specificity of expression within clusters. UMAP overlain with the normalized expression of TE markers (C) MER11A or (E) MER11C further confirms cluster specificity. (F) Mean expression for L1PA10 in the goblet cluster, and (G) UMAP of the goblet cluster with expression of L1PA10 overlain, separated by NIBD/CD. (H) Mean expression for SVA-D in the secretory progenitor cluster, and (I) UMAP of the secretory progenitor cluster with expression of SVA-D overlain, separated by NIBD/CD. P values were calculated by the Wilcoxon rank-sum test. ^#Adjusted P < .1, ∗adjusted P < .05.

Figure 7

CD disrupts stem cell homeostasis. (A) Differentially expressed gene within the stem cluster. Log₂ fold change is shown on the x-axis and the adjusted P value is shown on the y-axis. (B) Mean expression of 8 mediators of the Wnt pathway in NIBD and CD samples show a consistent decrease in response to CD. (C) Mean expression of CCDC115 and RNMT across clusters shows a significant decrease only occurs in the stem cluster. Percentage of cells positive for (D) LGR5 (left) and (E) SMOC2 (left) or the average normalized expression for (D) LGR5 (right) and (E) SMOC2 (right) in the stem cluster for NIBD or CD samples. Dots relate to individual sample value. (F) Uniform Manifold Approximation and Projection (UMAP) of the stem cluster with expression of LGR5 or SMOC2 overlain, separated by NIBD/CD. (G) UMAP of the stem cluster with LGR5 expression (red) co-overlain with either CD74 (top) or LONP1 (bottom) (both blue) expression. Right: Scale for expression is shown. (H) Percentage of cells positive for CD74 (left top) and LONP1 (left bottom) or the average normalized expression for CD74 (right top) and LONP1 (right bottom) in the stem cluster for NIBD or CD samples. Dots relate to individual sample value. (I) Shift in percentage of positive genes in either ISC-I, ISC-II, or ISC-III cell types between NIBD and CD. P values for expression bar plots were calculated using the Wilcoxon rank-sum test and the P value for the ISC subtype shifts was calculated using a Student t test. ∗P < .05, ∗∗P < .01, and ∗∗∗P < .001.

Figure 8

Subcluster analysis of SPIB+ cells shows distinct lineages. (A) Uniform Manifold Approximation and Projection (UMAP) after subclustering of SPIB+ cells with labeling of the 4 identified subclusters. (B) UMAP of SPIB+ cell subclusters overlain with gene expression from genes marking the entire cluster or distinct subclusters. (C) SPIB+ cell subcluster UMAP following cell type assignment using highly enriched markers of the subclusters. (D) Crypt–axis scores (low near crypt bottom, high near crypt top) of subcluster-assigned cells uncovers a colonocyte signature in the OTOP2+/BEST4+ subcluster. Clusters are arranged on the x-axis by mean crypt–axis score. The size of the dot corresponds to the level of expression of SELENOP, a known marker of the top of the crypt. (E) Mean cluster cell abundances across NIBD and CD samples shows a significant increase in OTOP2+/BEST4+ subcluster in CD. Dots show abundances of individual samples. (F) UMAP of either all clusters (top) or SPIB+ cell subclusters (bottom) overlain with the expression of markers of secretory progenitor cells. (G) Mean normalized expression of LYZ in NIBD and CD samples in the LYZ+ proliferating subcluster. (H) LYZ-detected immunofluorescence in colonic crypts from either NIBD (n = 181 crypts, 5 patients) or CD (n = 76 crypts, 2 patients) patient examples (left) and quantified by sample (right). Mean percentage of crypts positive for LYZ in NIBD or CD. (I) Lineage reconstruction as determined by partition-based graph abstraction (PAGA). P values are calculated by a Student t test. ∗P < .05. DAPI, 4′,6-diamidino-2-phenylindole.

Figure 9

Expression of CD-associated risk genes across clusters. (A) Dot plots showing expression and percentage of cells expressing CD-associated risk genes (columns) across clusters (rows). For each gene, the dot color corresponds to the average scaled expression and the size of the dot corresponds to the percentage of cells in the cluster expressing the gene. (B) Uniform Manifold Approximation and Projection (UMAP) of all epithelial clusters showing MUC2 (blue) and ITLN1 (red) expression. (C) UMAP of SPIB+ cell subclusters overlain with gene expression from NOTCH2 or ATG16L2. avg.exp.scaled, average expression, scaled; pct.exp, percent of cells expressing gene.

Figure 10

Differential expression of CD-associated risk genes across clusters. (A) Dot plots showing the log₂ fold change of CD compared with NIBD within clusters, the direction of change, and the statistical significance of CD-associated risk genes (columns) across clusters (rows). For each gene, the dot color corresponds to the log₂ fold change, the shape of the dot corresponds to the direction of change, and the color of the dot corresponds to the adjusted P value of the comparison. Adjusted P value was calculated from the Wilcoxon rank-sum test. DE, differential expression; padj, adjusted p-value.