Multimodal single-cell analyses reveal mechanisms of perianal fistula in diverse patients with Crohn's disease

Abstract

Background: Crohn's disease complicated by perianal fistulae is more prevalent and severe in patients of African ancestry.

Methods: We profiled single cells from diverse patients with Crohn's disease with perianal fistula from colorectal mucosa and fistulous tracts. Immunofluorescence was performed to validate predicted cell-cell interactions. Unstimulated monocytes were chronically cultured in diverse cohorts. A subset was analyzed by single-nucleus RNA + ATAC sequencing.

Findings: Fistulous tract cells from complete proctectomies demonstrated enrichment of myeloid cells compared to paired rectal tissues. Ligand-receptor analysis highlights myeloid-stromal cross-talk and cellular senescence, with cellular co-localization validated by immunofluorescence. Chitinase-3 like-protein-1 (CHI3L1) is a top upregulated gene in stromal cells from fistulae expressing both destructive and fibrotic gene signatures. Monocyte cultures from patients of African ancestry and controls demonstrated differences in CHI3L1 and oncostatin M (OSM) expression upon differentiation compared to individuals of European ancestry. Activating protein-1 footprints are present in ATAC-seq peaks in stress response genes, including CHI3L1 and OSM; genome-wide chromatin accessibility including JUN footprints was observed, consistent with reported mechanisms of inflammatory memory. Regulon analyses confirm known cell-specific transcription factor regulation and implicate novel ones in fibroblast subsets. All pseudo-bulked clusters demonstrate enrichment of genetic loci, establishing multicellular contributions. In the most significant African American Crohn's genetic locus, upstream of prostaglandin E receptor 4, lymphoid-predominant ATAC-seq peaks were observed, with predicted RORC footprints.

Conclusions: Population differences in myeloid-stromal cross-talk implicate fibrotic and destructive fibroblasts, senescence, epigenetic memory, and cell-specific enhancers in perianal fistula pathogenesis. The transcriptomic and epigenetic data provided here may guide optimization of promising mesenchymal stem cell therapies for perianal fistula.

Funding: This work was supported by grants U01DK062422, U24DK062429, and R01DK123758.

Keywords: AP-1; CHI3L1; Crohn's disease; PTGER4; Translation to patients; ancestry; myeloid; perianal fistula; senescence; single cell; stromal.

Figure 1.. Single cells from perianal fistula tracts implicate myeloid-stromal interactions.

A. Sampling strategy for acquiring biopsies from fistula tracts and paired inflamed distal rectal mucosa from proctectomies. B. H&E cross-sections of n=3 excised perianal fistula tracts, marked by dashed boxes. Patient ID in lower right corner. Scale bar = 5mm. C. Uniform manifold approximation and projection (UMAP) of annotated cells from n=3 paired fistula and rectum samples. D. Proportions of each cluster comprised by cells originating from fistula tract or rectal mucosal biopsies. N = cells in group. E. Heatmap of the total number of predicted ligand-receptor interactions between cell types (output in Table S2). Mono/Mac; Monocytes, Macrophages. F. Dot plot of significant ligand-receptor pairs predicted by CellPhoneDb, representing the likelihood that the mean value of a pair per cell-cell interaction is cell-type specific (see Table S2). Interactions on the x-axis are named in order of First_Partner|Second_Partner relative to the ligand/receptor pairs listed on the y-axis. G. Immunofluorescence staining showing proximity of ferritin (red) and SCARA5 (green) protein expression along the fistula tract lining (n=3 biological replicates) and corresponding sections of H&E-stained tissue. Dashed lines indicate the fistula tract lining. Scale bars for IF images: 200 μm; H&E images: 250 μm. H. UMAP of PerianalCD scRNAseq dataset by cluster annotation, n=15 patients. See also Figures S1–S2.

Figure 2.. Myeloid and stromal cell re-integration: gene expression from perianal fistula tracts.

A, B. UMAP of sub-clustered and re-integrated myeloid and stromal cells from colorectum and fistula from n = 15 PerianalCD cohort by cluster annotation (A) and by biopsy origin (B). C. Proportions of the total cells in each cluster originating from fistula or colorectal biopsies. Fibro, fibroblasts; FRC, fibroblastic reticular cells; MonoMacs, monocytes/macrophages. D. Projection of predicted ligand-receptor pairs in Fig 1F onto re-integrated myeloid-stromal clusters. E. Top DEGs in stromal cell clusters between fistula tract and colorectal biopsies. F. Enrichment of pre-defined fibrotic and destructive transcriptional modules in myeloid-stromal clusters. G. Top DEGs between integrated fibroblasts from fistula tracts and involved ileal biopsies, grouped by origin and expression of CHI3L1. H. Top DEGs in myeloid cell clusters between fistula tract and colorectal biopsies. I. Expression of macrophage activation path gene modules in myeloid cells by weighted kernel density estimation. Violin plots are annotated by the mean of the first module score component in fistula or tissue; Wilcoxon rank-sum test. J. Individual and joint weighted kernel density estimation of CD14 and CHI3L1 in myeloid and stromal clusters, as in Fig 2A. See also Figure S3.

Figure 3.. Fistulous tract and rectal immunofluorescence.

A, B. IF staining of fistula tract cross-sections with corresponding H&E sections for n = 3 fistula samples recruited for this study (A) and n = 2 fistula cross-sections from biobanked specimens (B); blue (cell nuclei), green (CHI3L1), red (CD14). Dashed lines indicate the fistula tract lining. Scale bars: 500 μm. C. Representative images of co-localized CHI3L1 and CD14 staining; scale bar = 20 μm. Images have been cropped. D. Quantitation of CHI3L1 and CD14 mean fluorescence intensity localized to the fistula tract lining vs. deeper tissue in the same cross-section. N=5 patients, each dot is 5-10 images quantified per section, 5 regions per deep vs. fistula location. Unpaired t-test. E, G. IF staining of CHI3L1 (green), IL13Rα2 (magenta) and α-SMA (red) in Perianal11 (E) and CD14 (red), CCL18 (green) and CD206 (magenta) in Perianal14 (G). Scale bars: 500 μm. Dashed lines indicate fistula tract lining. F. Quantitation of α-SMA and IL13Rα2 as in (D), n=3 patients. All mean fluorescence intensity values in arbitrary units. See also Figure S4.

Figure 4.. Distinct morphologic and molecular profiles of differentiated CD14+ monocytes and CHI3L1 from AA and EA individuals.

A. Experimental schematic of CD14+ cell 14-day in vitro culture. B. Representative images of day 14 cultures; scale bar 500 μm. Images have been cropped. C. Annotation of spindles in culture images by a trained AI model in a representative image. Scale bar: 500 μm. D. Quantitation of length, perimeter, area, and eccentricity measurements of day 14 cultured cells; one-way ANOVA with Tukey’s multiple comparisons test. N=8 controls, n=11 CD patients; 3 images per well of culture were captured, 3-10 wells per biological replicate. E. Cell marker genes measured by qRT-PCR relative to expression of housekeeping gene RPL32 on day 14 of culture; one-way ANOVA with Dunnett’s multiple comparisons test. Each point = 2-4 technical replicates, n=6 control, n=11 CD. F. Expression of CHI3L1 and PLAU relative to housekeeping on days 0 and 14 of culture; two-way ANOVA with Sidak’s multiple comparisons test. Each point = 2-4 technical replicates. G, H. Luminex assay of secretion of CHI3L1 (G) and OSM (H) by differentiated CD14+ cells, n=8 controls, n=11 CD patients, n=1 isolated fistula single cells in culture; n=2 technical replicates. Two-way ANOVA with Sidak’s and Tukey’s multiple comparisons tests. I. Expression of IL13RA2 and EGFR relative to housekeeping on days 0 and 14 of culture; two-way ANOVA with Tukey’s multiple comparisons test. Each point = 2 technical replicates. All error bars: standard error of the mean. * P ≤ 0.05 ** P ≤ 0.01 *** P ≤ 0.001 **** P ≤ 0.0001. See also Figure S5.

Figure 5.. Paired multiomic ATAC-sequencing and epigenetic regulation of fistula-associated regions.

A. Schematic for the parallel processing of colorectal biopsies for cytosolic scRNAseq and joint nuclear RNA + ATAC sequencing (10X Multiome). B. UMAP of annotated nuclei from the n = 6 Multiome cohort. C – F. Bias-corrected peaks in the promoters of CHI3L1 (C), OSM (D), PLAUR (E), and upstream enhancer of PLAU (F); each bin track is annotated with the number of footprints (bound transcription factor motifs) predicted. Grey arrow in (C) and grey bars in (D) indicate the genomic positions of the predicted bound footprints relative to the peaks. G. Genome-wide aggregate footprints and expression in binned scRNAseq clusters of JUN and FOS. For footprints, x-axis is the base pair (bp) position in either direction from the center of the motif (0), and the y-axis shows the signal of the mean Tn5 insertion rate. H. scVelo (RNA velocity) phase portrait of the ratio of spliced:unspliced FOS transcripts in scRNAseq clusters. Black dashed line indicates the estimated steady-state ratio. I. Significant regulons by z-score, calculated from the SCENIC AUCell score, from n=3 fistula scRNAseq samples. See also Figures S6, S7.

Figure 6.. Cell-type specificity of transcriptional regulation via multimodal integration.

A, B. Ranked region-based (A) and gene-based (B) eGRN specificity scores (eRSS) of n=26 high quality eGRNs in the Multiome bins. Activity of the eGRN (equivalent to eRegulon) is quantified by the eRSS on the y-axis. High quality is defined by correlation of >0.6 of both TF-gene and TF-region relationships. The +_+ symbols in each eGRN name correspond to the relationships between TF-gene and region-gene, respectively. C. The number of transcription factors within intersections of nucleus bins (vertical bars) and per set (individual bins, horizontal bars). Dots and connections signify the members of the intersect (vertical). Cyan represents TFs that are associated with IBD GWAS loci; gray bars are comprised by TFs with no known genetic association to IBD. See also Figure S7.

Figure 7.. Integration of top AA IBD genetic signals and chromatin footprints.

A. Total peaks, genome coverage, and overlap with significant IBD variants per nucleus bin. Fold change represents the difference in the % of the whole genome captured in ATAC peaks and the % of genome wide significant loci (1E-8) captured in peaks. B. Bias-corrected peaks in the region encompassing fine-mapped credible PTGER4 SNPs. Bins are annotated with the number of bound TFs (footprints) detected within the peak. Positions of the n=22 credible SNPs in this intergenic region are mapped below the plot (triangles). Above: motif plots for the two bound TF motifs closest to the lead SNP (rs11742570); RORC −7bp, RORA −9bp. C. Bias-corrected peaks in the PTGER4 gene and promoter. Bins are annotated with the number of footprints. D. Expression of PTGER4 in cytosolic scRNAseq from the PerianalCD cohort; UMAP clusters as in Figure 1H. E. Bias-corrected peaks in the region encompassing fine-mapped credible LACC1 SNPs; track annotations as in (B), SNP positions in the region in are indicated by the triangles below. F. Bias-corrected peaks in the LACC1 gene and promoter. Bins are annotated as above. G. Expression of LACC1 in cytosolic scRNAseq; UMAP clusters as in Figure 1H. See also Figure S7.