Mice
C57BL/6 J (B6) and Foxp3GFP, Rag1−/−, CD2DsRed, Foxp3CD2Il10GFP, UbPA-GFP and Nur77GFP mice were bred and maintained under specific pathogen free conditions in an accredited animal facility at the University of Oxford. Lta−/− mice were purchased from Jackson Laboratories. HH7-2tg TCR transgenic mice, referred to as TCRHh here, were provided by D. R. Littman. TCRHh mice were bred with CD2DsRed mice and either Foxp3CD2Il10GFP or Nur77GFP mice to generate TCRHhCD2DsRedIl10GFP or TCRHhCD2DsRedNur77GFP, respectively. TCRHh mice were bred to Rag1−/− to generate TCRHhRag1−/− mice. UBPA-GFP mice were bred to CD2DsRed and TCRHh mice to generate CD2DsRedUbPA-GFP and TCRHhCD2DsRedUbPA-GFP, respectively.
Mice were free of Helicobacter spp. and other known intestinal pathogens, were age- and sex-matched and between 6 and 12 weeks old. Animals were randomly assigned to experimental group, and cages contained mice of all different experimental groups. All experiments were conducted in accordance with the UK Scientific Procedures Act of 1986, and by persons holding a personal license. The project licence governing the mouse studies was reviewed by the University of Oxford’s Animal Welfare and Ethical Review Board and approved by the Home Office of His Majesty’s Government.
No statistical methods were used to predetermine sample size. Sample sizes were based on previous similarly designed experiments from our research group. The spatial transcriptomics experiment included four mice per group to balance statistical power with cost. For other experiments we aimed for a minimum of five mice per experimental group. Exact mouse numbers for each experiment are included in the figure legends. Mice were assigned to different experimental groups at random. Mice were co-housed and littermate when possible. Each cage contained all treatment conditions. Animal studies were not blinded. Histopathology scoring was conducted by two independent assessors, one of whom was blinded.
Flow cytometry
Mouse cells were stained with combinations of the following monoclonal antibodies, all purchased from Biolegend, Invitrogen, or eBioscience: CD4 (RM4-5), TCRβ (H57-957), CD45.1 (A20), CD45.2 (104), CD11c (N418), CD11b (M1/70), anti-human CD2 (TS1/8), CXCR5 (L138D7), PD-1 (J43), FOXP3 (FJK-16s), RORγt (Q31-378), Ki-67 (SolA15). Dead cells were excluded using efluor 780 fixable viability dye (eBioscience). For transcription factor staining, cells were stained with surface markers prior to fixation and permeabilization using the FOXP3 staining buffer kit (eBioscience) according to manufacturer instructions.
Immunofluorescence staining
Swiss-rolled caecum tissues were fixed overnight at 4 °C in PLP buffer (1% paraformaldehyde, l-lysine 0.2 M pH 7.4 and 32 mg NaIO4). Then, tissues were dehydrated in 20% sucrose for at least 4 h at 4 °C and embedded in OCT compound (Avantor). Seven-micrometre cryosections were rehydrated, blocked and permeabilized with PBS, 1% goat serum, 1% BSA, 0.3 M glycine, 0.3% Triton X-100 for 1 h at room temperature. Sections were stained with the following antibodies: Alexa Fluor 488 anti-mouse CD172a (SIRPα) (clone P84, 5 μg ml−1 Biolegend), Alexa Fluor647 anti-mouse CD11c (clone N418, 5 μg ml−1 Biolegend) or FITC/Alexa Fluor 594 anti-mouse CD206 (clone C068C2, 5 μg ml−1 Biolegend) and anti-mouse CD64 (clone X54-5/7.1, 4 μg ml−1 Biolegend). Sections were stained overnight at 4 °C. Before imaging, nuclei were counterstained with Hoechst. Images were acquired using Zen Blue software on a ZEISS 980 Airyscan inverted microscope equipped with a motorized stage. Diode laser lines were used for excitation: violet (405 nm), blue (488 nm), yellow (514 nm) and red (639 nm). All images were acquired with a 25× (NA 0.8) LD LCI Plan-Apochromat oil-immersion objective.
Colon tissue was embedded in OCT (Tissue-Tek) as Swiss rolls and sectioned at 7 μm. Slides were fixed with 3.7% formalin (Merck) and blocked with 10% donkey serum (Sigma Aldrich) and 1% Fc block (eBioscience) in permeabilization buffer (Foxp3/Transcription factor staining buffer set, eBioscience). B220 (RA3-6B2), CD4 (RM4-5), MHC Class II (M5/114.15.2), gp38 (8.1.1), IgD (11-26 c.2a) and BCL6 (IG191E/A8) (all Biolegend) were stained overnight in blocking serum.
Isolation of lymphocytes from spleen, lymph node and intestinal tissue
Intestinal tissues were washed twice in RPMI (Sigma Aldrich)/10%FCS/5 mM EDTA at 37 °C with agitation for 25 min to remove epithelial cells. CP and OLS were removed under 40× bright-field microscopy using a scalpel and a 16G needle and syringe, respectively. Remaining colon and caecum tissue, OLS and CP were digested for 40 min at 37 °C with agitation in RPMI/10% FCS/15 mM Hepes with 100 U ml−1 collagenase VIII (Sigma Aldrich) and 20 mg ml−1 DNase I (Sigma Aldrich). Leukocytes from colon and caecum tissue were recovered at the interface of a 40/70% Percoll gradient (Fisher Scientific).
Spleens and MLNs were mechanically disrupted, and splenic red blood cells were lysed with ACK lysis buffer.
Peripheral blood was collected by cardiac puncture and red cells were lysed with ACK lysis buffer.
Hh culture and oral gavage
Hh NCI-Frederick isolate 1 A (strain 51449) was grown on blood agar plates containing 7.5% laked horse blood (Thermo Scientific) and Skirrow Campylobacter supplement (Oxoid) under microaerophilic conditions at 37 °C with agitation. Cultures were expanded for 48 h in Tryptone Soy Broth (TSB, Fisher) containing 10% FCS (Gibco) and the above antibiotics. The concentration of bacteria was determined by optical density (OD) analysis at 600 nm. Mice were fed 1 × 108 colony-forming units of Hh (equivalent to 1 OD unit) by oral gavage using a curved 22G needle for a total of 2 doses 24 h apart.
Lymph node lymphocyte egress blocking experiment
Host mice were treated every 24 h with 1 mg kg−1 of the sphingosine-1-phosphate antagonist, Fingolimod (FTY720, Sigma Aldrich) via intraperitoneal injection at the indicated timepoints after naive TCRHh cell transfer.
Sorting and adoptive transfer of naive TCRHh T cells
Naive T cells were isolated from TCRHh mice splenocytes and sorted by flow cytometry as CD45+CD3+CD4+CD44lowCD62LhiVβ6+ (Extended Data Fig. 1d), with up to 2% contamination with nTreg cells. All monoclonal antibodies were purchased from Biolegend or eBioscience: CD3 (145-2C11), CD11b (M1/70), CD11c (N418), B220 (RA3-6B2), CD62L (MEL-14), TCRVβ6 (RR4-7), CD44 (IM7), CD45 (30-F11), CD4 (RM4-5). Sorted cells (2 × 103 or 5 × 104) were injected by intravenous injection into the tail vein for flow cytometric or in vivo live imaging respectively.
In vitro co-culture
Bone marrow stem cells were extracted from wild-type mouse femurs and cultured for 7 days in RPMI (Sigma) supplemented with 1% penicillin-streptomycin (Sigma), 10% FCS (Life Technologies), 1% Glutamax (Invitrogen) and 20 ng ml−1 mouse GM-CSF (Peprotech). Bone marrow-derived dendritic cells (BMDCs) were plated at a density of 1 × 104 cells per well overnight. Hh peptide (1 mg ml−1, Genscript) was added 1 h prior to plating 1 × 105 sorted naive TCRHhNur77GFP T cells in RPMI/10% FCS/1% Glutamax/1% penicillin-streptomycin and 50 mM β-mercaptoethanol (Life Technologies). Anti-mouse I-A/I-E antibody was added at 4, 12, 24, 48, 72 and 96 h after the plating of TCRHhNur77GFP T cells.
Two-photon microscopy
Mice were anaesthetized with isoflurane, the caecum exposed and immobilized with a suctioning imaging window46. Samples were illuminated with 910 nm <70 fs pulsed light from Mai-Tai laser and collected using a 20× water-dipping lens and the spectral detector of a Zeiss 880 multiphoton microscope (Carl Zeiss). Images were linearly unmixed using the Zen software (Carl Zeiss) to separate autofluorescence, collagen, eGFP, DsRed and Texas-red dextran based on single-colour controls. All in vivo live images of the LA and LP were performed in the caecum due to ease of access.
Intravital videos were drift-corrected based on mucus or collagen signal. Images were smoothed using a Gaussian filter for display.
Image analysis
Intravital microscopy images were analysed using Imaris 9 (Bitplane). Following unmixing, autofluorescence was subtracted from DsRed and GFP channels. GFP+ and DsRed+ cells were marked using the Surface Creation Wizard, and their co-expression and location within LP and LA compartments were recorded.
Immunofluorescence images were analysed using Imaris 10.0 (Bitplane). LAs were defined based on DAPI staining showing nuclear density and surrounding epithelial morphology, allowing for unbiased region selection. LP was chosen based on location and cellular density based on DAPI stain so that each region of interest collected and analysed for LA contained at least six LP surfaces of matching size. Surfaces for each cell type of interest were created using CD11c+ for BMDCs, CD206+ for CD206+ macrophages using the Surface Creation Wizard, which was applied to all images collected with the same conditions. CD11c+ surfaces were further subdivided based on median SIRPα staining levels within the surface.
Photo-activation
CD2DsRedUbPA-GFP hosts were colonized with Hh on day 0. Seven days later 50,000 naive TCRHhCD2DsRedUbPA-GFP cells were transferred into the Hh-colonized CD2DsRedUbPA-GFP hosts so both host and donor cells were photo-activatable. On day 21 the MLN and caecum were removed from 6 mice. The caecum was shaken at 37 C in RPMI + BSA + EDTA for 40 min and 20 min to remove the epithelium. The CP was separated from the rest of the tissue, and the remaining tissue was divided into one third without photo-activation as a negative control for FACS gating, one third for LA photo-activation, and one third for LP photo-activation.
Samples were maintained in RPMI + BSA+Hepes on ice in the dark for the duration of the experiment with cold media flowed over the tissue during photo-activation. A Zeiss 880 upright multiphoton microscope (Carl Zeiss, Germany) fitted with two tunable lasers (Mai-Tai tuneable BB laser 710–990 nm, pulse width <80 fs and Mai-Tai tuneable 690–1040 nm, pulse width <70 fs) and a 20× water-dipping objective was used for tissue photo-activation. The microscope was set to dynamically unmix GFP, DsRed, and collagen based on pre-collected single-colour controls. Samples were imaged with 910 nm light to identify regions of interest based on CD2DsRed. MLN and CP T cell zones were defined as the densest T cell regions without gaps to exclude B cell zones. LA regions were defined as a cluster of cells with a diameter of at least ten CD2DsRed cells. LP region was defined as a region containing CD2DsRed cells distal to LA. After ROI definition, the second laser was turned on at 740 nm while imaging live. GFP photo-activation was observed dynamically to ensure sufficient photo-activation without toxicity. For each mouse 3–10 regions were photo-activated for each tissue microniche. Each photo-activation region comprised approximately 40,000 µm3 of tissue (70 µm diameter × 10 µm depth).
After photo-activation, tissues were minced and digested for 30 min in RPMI+Hepes with 100 U ml−1 collagenase VIII (Sigma Aldrich) and 20 mg ml−1 DNase I (Sigma Aldrich). The digested tissue from the six mice was pooled per tissue microniche into one sample. Isolated cells were stained with efluor 780 fixable viability dye (eBioscience). Photo-activated GFP+ cells were sorted using a four laser BD FACSAria III based on gates defined by unactivated samples from the same mice (Extended Data Fig. 5b).
Single-cell RNA library construction and sequencing
For scRNA-seq experiments, the Chromium Single Cell 5′ version 2 reagent kit and Chromium Single Cell Mouse TCR Amplification Kit (10x Genomics) were used. Sorted cells were loaded onto each channel of the Chromium Chip K following the manufacturer’s instructions and the chip was inserted in the Chromium Controller for droplet encapsulation. cDNA synthesis, amplification, gene expression (GEX) and targeted TCR was performed on single cells, according to the manufacturer’s protocol (CG000331). Sequencing was performed on the Illumina Novaseq 6000 system. Gene expression libraries were sequenced at a targeted depth of 50,000 reads per cell, using the following parameters: Read1: 28 bp i7:10 bp, i5: 10 bp, Read2: 98 bp. TCR libraries were pooled at a ratio of 1:10 with the GEX libraries and sequenced at a target depth of 5,000 reads per cell.
scRNA-seq analysis
Pre-processing of 10x Genomics scRNA-seq and scTCR-seq data
scRNA-seq raw sequencing data were processed using the CellRanger “multi” software (version 6.1.1, 10x Genomics) with the mm10 2020-A mouse reference genome (official 10X mouse pre-built reference). Single-cell TCR-sequencing (scTCR-seq) data were aligned and quantified using the CellRanger ‘multi’ software (v.6.6.1) and the reference vdj_GRCm38_alts_ensembl-5.0.0 was used with default settings.
Quality control and processing of scRNA-seq data
Data pre-processing was performed using the ScanPy workflow (v. 1.8.1)47. ScanPy (v. 1.7.1), Anndata (v. 0.7.5), Pandas (v.1.2.3), NumPy (v.1.20.1), and Python (v.3) were used to pool single-cell counts and conduct downstream analysis. For each run, SoupX algorithm48 was run with default parameters to remove ambient mRNA from the count matrix. Doublet detection was performed using the Scrublet algorithm (https://github.com/AllonKleinLab/scrublet49) with percolation step, as previously described50. Additional doublet exclusion was performed throughout downstream processing based on unexpected co-expression of canonical markers, such as Cd3d (T cells) and Cd19 (B cells). Cells with fewer than 1,000 unique molecular identifier (UMI) counts, fewer than 600 detected genes and more than 15% mitochondrial reads were excluded from downstream analysis. Genes were filtered out for expression in less than three cells. Gene expression for each cell was normalized (scanpy.pp.normalize_total, scaling factor 10,000) and log-transformed (scanpy.pp.log1p). Downstream analyses included variable gene detection (scanpy.pp.highly_variable_genes) and data feature scaling (sc.pp.scale). Cell cycle score was calculated using the expression of the cell cycle genes in Supplementary Table 1. Cell cycle score, UMI counts, the percentage of mitochondrial reads, and the percentage of Ig reads (calculated based on the genes provided in Supplementary Table 1) were regressed out during scaling the data. Dimensionality reduction (scanpy.tl.pca, based on highly variable genes) and Leiden-graph-based clustering (scanpy.tl.leiden, with clustering resolution manually adjusted, 0.3–1.5) were carried out. Cell lineages were annotated on the basimarker gene expression for each cluster (sc.tl.rank_genes_groups, method = ‘wilcoxon’).
Cell-type annotation with CellTypist
CellTypist is a cell-type database, server and pipeline for automatic annotation of scRNA-seq developed at Teichman lab (https://github.com/Teichlab/celltypist, https://www.celltypist.org). To assemble a mouse intestinal reference dataset, scRNA-seq data were collected from eight publications covering different cell lineages from small and large intestine as well as one additional dataset of sorted B cells from spleen, to cover in detail the annotation of germinal centre B cell populations (Supplementary Data 2).
For each dataset, the raw count matrix was downloaded along with the accompanying cell meta information. After removing trivial cell types annotated by the original studies (for example, ‘doublets’ and ‘unresolved’), a total of 171,271 cells were obtained representing all major cell populations in the mouse gut. Cell-type names from different datasets were next standardized to achieve a common scheme of nomenclature. Specifically, the similarity of transcriptomes between each pair of cell types across datasets was examined and the two cell types were merged only if they corresponded to a single cell type (for example, ‘enteroendocrine cell’ from the Tabula Muris was renamed to ‘enteroendocrine’ as was used in Biton et al.30). Finally, after in-depth inspection, 126 cell types were harmonized from the eight datasets. A CellTypist model then was created based on logistic regression classifiers, as described in detailed48. The model is publicly available at https://celltypist.cog.sanger.ac.uk/models/Mouse_Gut_Casado/v2/Adult_Mouse_Gut.pkl. Cell identities were predicted using the resulting model, followed by manual curation. Cells from each lineage were further subclustered and Leiden clustering was repeated for fine-grained annotation of the cell types and states. A description of cell-type annotations for each lineage is provided in Extended Data Fig. 5c and Supplementary Data 3 and 4. The differentially expressed genes for the cell types in each lineage can be found in the Supplementary Table 1.
For prediction on cycling regulatory T cells (prolif. Treg cells), eTreg and cTreg cells were used as a training reference. The model was built applying default parameters, and prediction was performed without majority voting.
scTCR-seq downstream analysis
The Python package scirpy (v. 0.12.2)51 was used to extract the V(D)J sequence information from the CellRanger output file filtered_contig_annotations.csv. Productive TCRαβ chains were determined using the scirpy.tl.chain_qc function, and cells without V(D)J data or with two pairs of productive TCRαβ chains were removed from the analysis. Clonotypes were defined with the function scirpy.tl.define_clonotypes based on the CDR3 nucleotide sequence identity and the V-gene usage for any of the TCR chains (either VJ or V(D)J need to match). For cells with dual TCRα or TCRβ chains, any of the primary or secondary receptor chain matching was considered for the clonotype definition. Clonotype networks were constructed using the pp.ir_dist function to compute distances between CDR3 nucleotide sequences (using identity as a metric) and tl.define_clonotype_clusters function to designate the clonotype clusters, removing the clonotypes with less than two cells. The TCR metadata were combined with the transcriptome data for downstream analysis and comparison of different T cell populations. Hh-specific TCR data were retrieved from Xu et al.15, and TCR sequences were obtained from https://www.ncbi.nlm.nih.gov/nuccore and mapped using the IMGT/V-QUEST alignment tool52. Hh.7-2 transgenic TCR (tgTCR) clonotypes were identified by expression of TRAV9-1/TRBV19 gene segments and the TCRβ CDR3 amino acid sequence, including those clonotypes with missing TCRα chain.
RNA velocity analysis
RNA velocity analysis53 was performed using the scVelo (v.0.2.4) package [10]. RNA velocity was estimated by distinguishing unspliced and spliced mRNAs using the velocyto package (v.0.17) (https://velocyto.org/velocyto.py/54). Data were subclustered on Treg cell, filtering out the subsets with fewer than ten cells per gut region (that is, excluding eTregs_MLN, eTregs_CP, cTregs_LA, cTregs_LP and Prolif-Tregs_CP). The dataset was then merged with the velocyto output (merged loom files) and pre-processed for detection of minimum number of counts, filtering and normalization (scvelo.pp.filter_and_normalize). The functions scvelo.pp.moments, scvelo.tl.velocity and scvelo.tl.velocity_graph were used to compute velocities using the stochastic mode in scVelo. The function scvelo.pl.velocity_embedding_stream was used to project the velocity information onto the UMAP. To test which genes have cluster-specific differential velocity expression and visualize selected genes, the functions scvelo.tl.rank_velocity_genes and scvelo.pl.velocity were applied. Velocity pseudotime was calculate with the function scvelo.tl.velocity_pseudotime.
Cell-type scoring
A list of mouse genes involved in the MHCII complex (Supplementary Table 2) was used for surface MHCII scoring. Cells were scored using the scanpy.tl.score_genes function according to the expression values of all genes.
Cell–cell communication analysis
The CellPhoneDB55,56 Python package (v.3 .0) was used to infer putative cell–cell interactions. The scRNA-seq dataset was split by gut region and cell types with <30 cells in a given region were filtered out. Human-mouse orthologue genes were retrieved using the biomaRt package57, and only one-to-one orthologous genes were considered. CellPhoneDB was applied on the normalized raw counts and fine cell-type annotations of myeloid, T cells and ILCs from LA and LP (separately for each gut region), using default parameters. To identify the most relevant interactions, specific interactions of Treg cells with myeloid cells and T cells/ILCs were selected and filtered for the ligand–receptor pairs that were significant (P ≤ 0.01) and ‘curated’. The selected interactions were plotted as expression of both ligands and receptors in relevant cell types. The ktplots R package (https://github.com/zktuong/ktplots/tree/plot_cpdb3; https://doi.org/10.5281/zenodo.5717923) was used to visualize the significant interactions per cell-type pair using a chord diagram.
Fresh frozen Visium sample preparation
Caecum and proximal colon tissue from Hh-infected and Hh/anti-IL10R treated mice were removed, cut longitudinally and cleaned of stool with cold phosphate buffered saline (PBS). The tissue was positioned luminal side up, and rolled into a Swiss roll from the caecum to the proximal colon. The tissue was placed into a histology plastic cassette and snap frozen for 1 min in dry-ice-cooled isopentane. The frozen tissue was embedded in OCT on dry ice and stored at −80 °C. The samples were selected based on tissue morphology and orientation (H&E-stained sections) and RNA integrity number, obtained using High sensitivity RNA ScreenTape system, (Agilent 2200 TapeStation). OCT blocks were sectioned at 10 μm thickness in a −20 °C cryostat (Leica CX3050S) at 10 μm, and transferred onto a 6 mm2 capture area on a Visium 10X Genomics slide. Visium spatial tissue optimization was performed, and an optimal permeabilization time of 24 min was selected. The Visium slides were processed according to manufacturer’s instructions, before fixing and staining H&E for imaging. H&E-stained slides were imaged at 40× on Hamamatsu NanoZoomer S60. After transcript capture, sequencing libraries were prepared according to the 10X Genomics Visium Spatial Transcriptomic protocol and sequenced on the Illumina Novaseq 600 system.
Visium spatial transcriptomics data analysis
10x Genomics Visium sequencing data processing
After sequencing, 10x Genomics Visium spatial samples were aligned to the mouse transcriptome mm10 2020-A reference (as the scRNA-seq samples) using 10x Genomics SpaceRanger version 2.0.0. and exonic reads were used to produce mRNA count matrices for each sample. SpaceRanger was also used to align paired histology images with mRNA capture spot positions on the Visium slides. A custom image-processing pipeline was used for alignment of Visium slides and identification of the spots contained in the tissue, as described in ref. 58. Spots with fewer than 500 UMI counts, and more than 15% mitochondrial genes were removed from the analysis. Data from different samples were concatenated and SCVI was used for batch correction59.
Spatial mapping of cell types using cell2location
To spatially map intestinal cell types defined by single-cell transcriptomics in the Visium data we used cell2location41. First, to obtain a complete single-cell reference of cell types and cell states in the mouse intestine we integrated our NICHE-seq data with 3 publicly available datasets of intestinal epithelial cells30, immune cells15 and enteric nervous system27 (Fig. 5a). Redundant cell annotations across different datasets were harmonized and curated manually. This scRNA-seq reference (untransformed and unnormalized mRNA counts) was then used in the cell2location pipeline, as described in detail previously43. In brief, reference signatures of cell states (63 cell populations) were estimated using a negative binomial regression model provided in the cell2ocation package. The inferred reference cell-state signatures were used for cell2location mapping that estimates the abundance of each cell state in each Visium spot by decomposing spot mRNA counts. The cell2location spatial mapping was done separately for Hh and Hh/anti-IL10R sections. The paired H&E images were used to determine the average number of cells per spot in the tissue (set to 30) and used as a hyperparameter in the cell2location pipeline. Cell-state proportions in each Visium spot were calculated based on the estimated cell-state abundances.
Two methods were used to identify cellular microenvironments in the tissue: manual annotation and conventional NMF analysis. Regions for manual annotation were defined based on H and E images. LA were defined by cellular density, whereas LP regions included histologically distinct areas both proximal and distal to the LA. NMF, implemented in the cell2location pipeline, was performed on cell abundance results by cell2location on each condition separately (Hh and Hh/anti-IL10R). The NMF model was trained for a range of factors and tissue zones (number of factors: n_fact) N = {5,…,30} and the decomposition into 18 factors was selected as a balance between segmenting relevant tissue zones (muscle compartment, lymphoid structures, lymphatics) and over-separating known zones into several distinct factors (Extended Data Fig. 8c).
Cell-state spatial enrichment analysis
Spots containing lymphoid aggregates and adjacent LP were manually annotated using the paired histology images of the spatial data in the 10x Genomics Loupe software. Cell-state proportions per spot were calculated based on the estimated abundances from cell2location and cell-state enrichments (odds ratio) in each manually annotated region were calculated as described60. In brief, the odds of target cell-state proportions were divided by the odds of the other cell-state proportions. Odds of cell proportions were calculated as the ratio of cell proportion in the spots of a structure of interest to that in the other spots. Statistical significance was obtained by chi-square analysis (scipy.stats.chi2_contingency) and the P value was corrected with the Benjamini–Hochberg method.
Spatial co-occurrence analysis of cell types
We quantified the degree of co-occurrence between cell types on the basis of their proportions inferred by cell2location. Specifically, since each cell type had an estimated abundance distribution across spatial spots, we calculated the Pearson correlation coefficient for any two given cell types as their co-occurrence rate. This calculation was conducted for each sample separately. Next, the log2 ratio between four Hh-only (control) samples and four Hh/anti-IL10R samples was defined as their fold change in co-occurrence rate and the significance (that is, P value) was assessed by a two-sided Student’s t-test.
Spatial ligand–receptor analysis
Ligand–receptor analysis on Visium data was performed using the Cell2location cell-type abundances and the stLearn package61 (https://github.com/BiomedicalMachineLearning/stLearn). In short, the connectomeDB2020_lit database for mouse was used as a reference of candidate ligand–receptor pairs. The st.tl.cci.run function was used to calculate the significant spots of ligand–receptor interactions within spot mode (distance = None), filtering out any ligand–receptor pairs with no scores for less than 20 spots, and using 10,000 random pairs (n_pairs). P values were corrected with the st.tl.cci.adj_pvals function using false discovery rate, Benjamini–Hochberg (adj_method = ‘fdr_bh’) adjusting by the number of spots tested per ligand–receptor pair (correct_axis = ‘spot’). Spot P values were displayed for particular ligand–receptor pairs (‘Ccl8_Ccr2’ and ‘Cxcl9_Cxcr3’) in the spatial context using the function st.pl.lr_result_plot.
Statistical analysis
Statistical analysis was performed using Prism 8 (GraphPad). t-Tests were used to compare two unpaired samples. For more than two groups, the ordinary one-way ANOVA was used. No samples were excluded from analysis. Mean with standard deviation shown unless otherwise indicated. Differences were considered statistically significant when P ≤ 0.05. Significance is indicated as *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001 and ****P ≤ 0.0001. Technical replicates were processed and analyses on the same day. Biological replicates are from independent experiments.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.